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SECTION  1 
SUMMARY 


The  a  priori  prediction  of  multinozzle  rocket  exhaust  flow  fields  is 
addressed  in  detail.  The  fundamental  requirements  for  accurate  prediction  of 
plume  IR  signature  are  derived  and  new  quantitative  relationships  between 
optical  siqnal  and  plume  properties  are  derived.  It  is  shown  that,  amonq  a 
variety  of  requirements,  plume  models  must  include  an  accurate  detailed 
description  of  the  three  dimensional  near  field  of  the  multinozzle  plume  self 
impingement  to  achieve  the  accuracy  and  reliabilty  of  the  optical  predictions, 
over  the  desired  altitude  range.  The  qualitative  structure  of  these  complex 
three  dimensional  flow  fields  is  explained  for  the  first  time.  Several  of  the 
regulating  flow  processes  thus  identified  are  three  dimensional  in  nature  and 
have  no  counterparts  in  classical  two  dimensional  supersonic  flow  theory.  One 
such  process,  the  intersection  of  two  three  dimensional  shock  surfaces,  is 
discussed  in  detail  and  a  qualitative  account  of  the  developing  pattern  is 
given.  A  numerical  procedure,  "the  floating  fitted  shock"  technique,  fit  the 
requirements  of  accuracy  and  generality  necessary  for  the  computation  of  the 
multinozzle  plume  flow  fields.  This  method  is  conceded  to  be  the  most 
desirable,  albeit  most  complex,  for  the  solution  of  supersonic  flows.  A 
computer  code  was  devised  which  contained  discrete  discontinuities  including 
slip  surfaces,  a  shock  surface  and  a  complex  sonic  shock/centered  expansion 
singular  point  based  on  a  boundary  point  calculation  which  properly  accounts 
for  the  three  dimensional  propagation  of  characteristics.  The  code  was 
successful  for  simplified  geometries  but  could  not  be  increased  in  generality 
to  handle  the  complete  flow  pattern.  An  analysis  for  the  Mach  disc  flow  field 
in  an  axisymmetric  plume  was  derived  which  leads  to  a  basic  interaction 
equation.  A  numerical  procedure  for  solving  this  equation  along  with  the 
other  governing  one  dimensional  equation  uses  the  novel  approach  of  first 
locating  the  sonic  throat  position  and  then  integrating  the  equations 
upstream.  Several  test  cases  are  presented  which  include  viscous  mixing  for 
the  first  time.  This  method  should  provide  more  reliable  calculation 
procedure  for  these  flow  fields  than  now  exists. 
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INTRODUCTION 

Detailed  understanding  of  jet  and  rocket  engine  exhaust  flow  fields  is 
required  in  a  wide  variety  of  Air  Force  space  and  missile  programs. 

Predictions  of  infrared  signature,  radar  cross  section,  electromagnetic  wave 
attenuation,  and  production  and  dispersion  of  noxious  pollutants  are  examples 
of  the  exhaust  system  properties  which  are  fundamental  in  both  conceptual 
systems  studies  and  actual  design  and  development  programs.  These  exhaust 
properties  are  the  subject  of  a  broad  area  of  study  known  as  plume 
phenomology-  a  multidisciplinary  study  encompassing  the  sciences  of  fluid 
mechanics,  chemical  kinetics  and  optical  radiative  transport  theory.  Often 
system  design  requirements  and  programmatic  study  definitions  require  state  of 
the  art  or  perhaps  beyond  predictive  capabilities  in  each  of  these  disciplines 
to  provide  the  desired  information  and  definitions.  This  research  program  was 
aimed  at  extending  the  state  of  the  art  of  fluid  mechanical  prediction 
techniques  while  keeping  in  mind  how  these  advances  would  fit  in  the  broader 
overall  plume  phenomology  program. 

Plume  fluid  mechanics  is  the  underlying  physical  science  in  plume 
phenomology  as  the  spatial  distribution  of  thermochemical  properties  it 
defines  provides  the  structure  upon  which  chemical  kinetics  and  subsequently 
radiative  transport  are  predicted.  The  main  area  of  interest  was  the  fluid 
dynamics  of  multiple  nozzle  exhaust  systems.  This  study  was  aimed  at  the 
inviscid  structure  of  this  flow  as  it  forms  the  "skeleton"  upon  which  is  built 
the  complete  flow  field  including  the  turbulent  (reacting)  mixing  layer  and 
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viscous  far  field.  This  portion  of  the  effort  was  therefore  the  study  of 
inviscid  three  dimensional  supersonic  flow  fields  containing  complex  shock 
systems.  Plume  flow  fields  can  contain  Mach  discs  and  hence  regions  of 
subsonic  flow.  Understanding  and  modelling  of  these  phenomena  was  a  second 
area  of  study. 

The  flow  field  created  by  the  exhaust  of  a  multiple  nozzle  exhaust  system 
is  a  complex  three  dimensional  flow.  The  inviscid  flow  defines  the  shock  wave 
structure  of  a  plume  which  is  of  prime  importance  in  predicting  plume 
observables.  The  shock  waves  are  responsible  for  both  local  sharp  increases 
in  temperature  and  pressure  and  far  field  temperature  increments.  The  far 
field  effect  is  a  product  of  the  entropy  rise  (total  pressure  loss), 
associated  with  the  shock  waves,  which  persists  downstream  showing  up  there  as 
an  increment  in  temperature  above  the  isentropic  far  field  temperature.  Both 
optical  radiation  and  chemical  kinetic  processes  are  governed  by  eauations  all 
of  which  contain  "Arrhenius"  type  exponential  factors  [exp  (-8/T)].  In  the 
case  of  optical  radiation  B  is  the  second  radiation  constant  divided  by  the 
wavelength.  The  "characteristic"  temperature  B  of  these  processes  is 
generally  high  so  that  chemical  activity  and  radiative  source  terms  are  most 
prominant  in  regions  of  high  temperature.  In  cases  where  B  is  much  larger 
than  the  maximum  temperature  the  regions  of  high  chemical  reaction  or 
radiative  emission  reduce  to  extremely  thin  sheets.  Thus  the  shock  wave 
structure  which  is  the  primary  factor  in  determining  both  local  and  far  field 
temperature  levels  is  a  central  determining  factor  in  chemical  activity  and 
optical  radiation.  For  example,  it  is  clear  from  many  axisymmetric  flow  field 
studies  (c.f.  Ref.  2-1)  that  shock  structure  and  peaks  in  plume  IR  station 
radiation  are  highly  correlated.  In  Ref.  2-2  the  central  role  of  far  field 
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temperature  on  plume  radiation  was  demonstrated  clearly.  It  is  for  these 
reasons  that  detailed  predictive  capabilities  of  the  three  dimensional  flow 
and  shock  structure  were  sought  for  the  multiple  nozzle  flows. 

The  pursuit  of  solutions  to  the  three  dimensional  plume  flow  fields 
requires  the  numerical  solution  of  the  Euler  Equations.  For  supersonic  flows 
these  are  a  set  of  hyperbolic  partial  differential  equations  which  forms  an 
initial  value  problem.  These  problems  are  suited  (in  theory)  to  straight 
forward  marching  numerical  solution  techniques.  Experimental  evidence  (Ref. 
2-3)  showed  that  there  were  at  least  three  and  most  probably  more  shock  wave 
configurations  possible  for  twin  nozzle  flow  fields.  The  specific 
configuation  being  determined  by  individual  nozzle  exit  properties,  nozzle 
spacing  and  background  pressure  and  flow  properties.  The  full  diversity  of 
flow  configurations  possible  may  be  very  broad  and  include  many  configurations 
not  yet  observed.  For  two  reasons  it  was  deemed  desirable  to  seek  numerical 
solutions  of  these  flow  fields  in  which  the  shock  waves  were  "fitted".  That 
is,  in  the  calculation  procedure  the  shock  wave  surfaces  are  considered  to  be 
discrete  discontinuous  surfaces  which  are  tracked  as  part  of  the  solution. 

The  first  reason  derives  from  the  discussion  in  the  previous  paragraph. 

Optical  and  radiative  transport  equations  are  extemely  sensitive  to 
temperature  so  that  the  most  accurate  temperature  predictions  are  desireable. 
Secondly,  the  alternative  approach  "shock  capturing"  produces  a  shock  wave  in 
a  two  dimensional  flow  that  is  portrayed  in  the  solution  as  a  dispersed 
compression  spread  over  three  to  five  mesh  intervals.  There  are  generally 
overshoots  and  undershoots  involved  in  the  flow  quantities  and  the  solution  in 
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the  neighborhood  of  the  shock  wave  is  therefore  unreliable  and  contains  order 
unity  inaccuracy.  The  multi nozzle  plume  is  a  three  dimensional  flow 
containing  shock  intersections  which  would  likewise  have  five  by  five  mesh 
regions  of  questionable  accuracy.  Therefore  it  appears  that  with  the  complex 
shock  structure  that  exists  in  the  subject  flow  "shock  capturing"  would 
produce  results  that  either  (1)  had  an  inordinate  percentage  of  mesh  points 
inaccurate  because  of  nearby  shocks  or  shock  intersections  or  (2)  required  an 
inordinate  number  of  mesh  points  to  cirmcumvent  this  problem. 

The  only  previous  complex  three  dimensional  calculations  (Ref.  2-4) 
employing  fitted  shock  waves  was  for  space  shuttle  type  configurations  in 
which  the  general  shape  and  topology  of  the  shock  pattern  is  known.  In  that 
case  very  elegant  precise  mapping  techniques  were  brought  to  bear  on  the 
problem  to  simplify  the  computer  program  logic.  Since  the  shock  wave  surface 
configurations  for  the  multinozzle  plume  flow  fields  is  not  known  precisely 
and  can  take  on  any  one  of  a  number  of  general  configurations  a  more  general 
numerical  technique  was  pursued.  This  is  referred  to  as  the  floating  "fitted" 
shock  wave  technique  and,  in  fact,  employed  floating  discontinuties  and 
singularities  more  general  than  shock  waves.  In  this  method  (cf.  Moretti 
Ref.  2-5)  the  shock  wave  surfaces  are  not  mapped  to  the  boundaries  of 
computational  domains,  rather  they  are  permitted  to  traverse  a  relatively 
stationary  computational  grid.  In  theory  the  scheme  does  not  require  a  priori 
knowledge  of  the  shock  wave  configuration  and  so  would  be  ideal  for  the 
problem  at  hand.  The  price  for  this  generality  is  very  heavy  and  two  fold. 
First  the  computer  logic  is  extremely  complex  as  it  must  anticipate  a  large 
number  of  geometrical  configurations  and  combination  of  configurations.  And, 


secondly,  the  details  of  the  three  dimensional  phenomenon  must  be  understood 
so  that  local  solutions  can  be  incorporated  into  the  flow  field.  For  example, 
in  this  "fitted  singularity"  approach  the  exact  local  solution  for  the  shock 
intersection  with  the  plume  boundary  was  incorporated  in  the  solution.  The 
exact  details  of  the  reflected  Prandtl  Meyer  expansion  wave  and  sonic  nature 
of  the  impinging  shock  are  employed  in  a  special  cell  calculation.  This 
singular  point  was  free  to  traverse  the  computational  mesh  as  the  solution 
progressed  and  pass  from  one  cell  to  the  next.  In  this  way  the  calculation 
procedure  is  similar  to  finite  element  methods  with  the  added  feature  that 
the  cell  containing  this  singularity  changes  automatically  as  the  shock  moves 
and  the  calculation  proceeds.  Unfortunately  there  are  several  other 
singularities  which  are  not  presently  understood  that  must  be  modelled  to 
accurately  predict  the  three  dimensional  flow  fields.  One  of  these  occurs 
when  two  three  dimensional  shock  surfaces  intersect  and  subsequently  develop 
into  an  irregular  reflection.  This  problem  was  studied  theoretically  and  a 
qualitative  description  of  the  flow  field  development  was  derived  employing 
hodograph  techiques. 

The  solution  of  plume  flow  fields  with  Mach  discs  open  up  an  entirely  new 
set  of  requirements  for  numerical  calculations  and  theoretical  assessment. 

The  exhaust  gas  flow  fields,  aside  from  regions  downstream  of  Mach  discs,  are 
supersonic  and  as  such  are  governed  by  hyperbolic  partial  differential 
equations.  These  equations  have  solutions  which  at  any  point  can  be  expressed 
solely  on  the  basis  of  flow  properties  upstream  of  that  point.  Therefore,  the 
numerical  solution  of  these  equations  proceeds,  at  least  in  principle,  by  a 
step  by  step  or  marching  procedure.  The  flow  downstream  of  Mach  discs  on  the 


other  hand  is  subsonic  and  as  such  is  governed  by  elliptic  partial 
differential  equations.  The  solution  to  these  equations  at  any  point  depends 
on  the  solution  at  all  neighboring  points  both  upstream  and  downstream.  Thus 
the  solution  at  any  specific  point  is  related  to  all  the  boundary  values  of 
the  flow  properties  surrounding  the  subsonic  region.  The  plume  containing  a 
Mach  disc  is  thus  a  mixed  type  flow  containing  regions  of  both  subsonic  and 
supersonic  flow  that  is  similar  in  many  ways  to  more  familiar  transonic  flows. 
These  flows  are  computed  numerically  by  overall  relaxation  schemes  (cf.  Ref. 
2-6.). 

The  Mach  disc  flow  field  has  a  critical  property  that  was  exploited  by 
Abbett  (Ref.  2-7)  to  explain  the  determining  factor  for  the  shock  triple  point 
location.  He  observed  that  the  flow  which  passes  through  the  normal  shock 
portion  of  the  Mach  disc  is  subsequently  accelerated  downstream  to  supersonic 
velocities.  This  subsonic/supersonic  stream  tube  is  analogous  to  a  De  Laval 
nozzle  with  a  choking  or  saddle  point  singularity  condition  at  the  throat. 

The  location  of  the  Mach  disc  itself  is  determined  as  that  position  which  is 
compatible  with  a  smooth  acceleration  through  the  sonic  singularity.  An 
important  approximation  which  greatly  simplifies  the  solution  of  plumes  with 
Mach  discs  was  given  by  Abbett  (Ref.  2-7)  and  later  employed  by  Salas  (Ref 
2-8).  The  subsonic  portion  of  the  flow  field  is  approximated  using  one 
dimensional  flow  analysis.  The  location  of  the  Mach  disc  is  estimated  and  the 
subsequent  supersonic  outer  flow  and  subsonic  one  dimensional  flow  is 
calculated.  The  unique  solution  is  determined  by  finding  the  location  of  the 
Mach  disc  which  results  in  sonic  velocity  occurring  at  the  same  point  as  the 
minimum  area.  Several  possible  Mach  disc  locations  are  employed  to  iterate 
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and  determine  a  bracket  for  the  possible  location  of  the  Mach  disc.  This 
method  has  been  used  (Salas  Ref. 2-8)  with  success  in  the  past  to  determine  a 
range  of  plume  flow  fields  with  Mach  discs.  However,  the  method  has  several 
drawbacks.  It  does  not  provide  any  theoretical  explanation  of  the 
relationship  of  the  solution  to  overall  flow  properties,  it  does  not  include 
viscous  effects,  and  it  has  questionable  numerical  reliability. 

The  forward  integration  of  the  subsonic  equations  is  a  tricky  procedure 
at  best  and  is  questionable.  The  set  of  equations  governing  the  subsonic 
streamtube  has  a  positive  eigenvalue.  Thus  it  possesses  exponentially  growing 
solutions  (unstable).  Compounding  this  fact,  the  exponential  factor  is 
proportional  to  (1  -  M^)“l  so  that  as  the  sonic  singularity  is  approached 
forward  marching  numerical  procedures  become  useless.  In  view  of  the  rapid 
changes  in  the  solution  near  the  throat  a  compromise  must  be  struck  between 
accuracy  and  the  ability  to  generate  solutions  at  all.  In  the  present  work  a 
solution  procedure  is  developed  which  resolves  this  dilemma.  In  the  inviscid 
case  the  sonic  throat  location  is  determined  first  and  then  the  equations  are 
solved  by  integration  in  the  upstream  direction. 

The  analysis  which  is  developed  leads  to  a  central  qoverninq  interaction 
equation.  This  equation  coupled  with  the  familiar  one  dimensional  flow 
equations  and  the  supersonic  flow  equations  for  the  outer  plume  stream  help  to 
delineate  the  underlying  processes  driving  the  solution.  The  stronq 
interaction  between  the  subsonic  flow  which  passed  through  the  normal  shock 
portion  of  the  Mach  disc  and  the  supersonic  flow  which  passed  through  the 
oblique  shock  is  evident  through  a  term  which  is  related  to  the  Reimann 
invariant  on  the  downward  running  characteristics.  Thus  the  explicit  effect 
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of  the  plume  outer  boundary  properties  on  the  Mach  disc  solution  is  made  clear 
for  the  first  time.  In  addition  viscous  mixing  between  the  two  streams  is 
incorporated  into  the  analysis  via  a  simple  one  dimensional  analysis.  In  the 
inviscid  case  determination  of  the  exact  location  of  the  sonic  throat  is 
possible  based  on  an  estimated  supersonic  flow  and  an  assumed  location  of  the 
shock  triple  point.  In  the  viscous  case  a  further  iteration  is  necessary. 

In  the  following  section  a  qualitative  analysis  of  the  multinozzle  plume 
flow  field  is  discussed.  The  various  three  dimensional  features  of  the  flow 
are  pointed  out  in  the  context  of  the  overall  flow  structure.  Section  4 
reviews  the  relationship  between  the  plume  fl uid/ thermochemical  distributions 
and  the  resulting  IR  signature.  The  mathematical  relationship  between  optical 
emission  and  flow  properties  is  derived.  In  Section  5  the  floating  fitted 
shock  numerical  technique  is  outlined  and  the  three  dimensional  boundary  point 
calculation  is  described.  A  sample  calculation  for  a  simplified  geometry  is 
presented.  Section  6  presents  the  study  of  the  intersection  of  two  three 
dimensional  shock  surfaces.  Section  7  briefly  describes  an  alternate  new 
method  for  finite  difference  calculations  which  might  reduce  the  enormous 
logic  load  on  the  fitted  discontinuity  programs.  Section  8  contains  the 
discussion,  analysis  and  numerical  computations  for  the  Mach  disc  flow. 
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SECTION  3 


STRUCTURE  OF  MULTINOZZLE  PLUME  FLOW  FIELDS 

The  structure  of  the  multinozzle  plume  flow  is  dominated  by  complex  three 
dimensional  flow  phenomenon.  These  processes  are  not  simply  the  extension  of 
familiar  two  dimensional  supersonic  flow  situations  into  a  third  dimension. 
Rather  they  are  new  and  peculiar  to  three  dimensional  flows  and  as  such  are 
basically  unknown  to  analysts.  Therefore  the  ability  to  accurately  model  and 
devise  numerical  schemes  and  their  associated  computer  codes  rests  heavily  on 
first  developing  some  understanding  in  these  areas.  There  are  two  basic 
situations  that  can  be  identified.  The  first  and  most  striking  is  the  problem 
of  the  intersection  cf  three  dimensional  shock  wave  surfaces.  The  multinozzle 
plume  flow  field  contains  several  shock  wave  surfaces  and  these  invariably 
intersect.  The  subsequent  development  of  the  shock  pattern  is  complex  and 
will  be  discussed  in  detail  in  a  later  section.  Another  situation  peculiar  to 
three  dimensional  flows  occurs  when  there  is  an  abrupt  change  in  geometry  or 
topology  in  the  flow  (i.e.  transition  of  a  shock  reflection  from  a  regular 
reflection  to  a  Mach  reflection).  This  brings  about  an  initialization  problem 
which  is  analogous  to  the  flow  at  the  leading  edge  of  a  wedge  or  the  point  of 
a  cone.  These  latter  two  are  well  known  two  dimensional  supersonic  flow 
situations  which  have  cataloged  solutions.  We  call  on  our  knowledge  of  these 
catalogs  to  initialize  or  reinitialize  two  dimensional  (or  axi symmetric )  flow 
fields  when  there  is  an  abrupt  change  in  geometry  (wall  angle).  These 
catalogs  do  not  exist  in  three  dimensional  flows.  Beyond  that  the  nature  of 
the  solution  is  not  known  in  any  of  these  situations. 
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In  order  to  study  the  structure  of  the  multi  nozzle  plume  flow  field  the 
structure  of  the  individual  undisturbed  exhaust  plumes  is  first  examined.  Each 
exhaust  plume  is  axi symmetric  and  underexpanded.  Figure  3-1  is  a  schematic  of 
one  of  these  plumes.  Supersonic  exhaust  flow  leaves  the  rocket  nozzle  at  the 
exhaust  plane.  The  ambient  pressure  at  the  exit  plane  is  lower  than  the 
exhaust  plane  pressure  and  the  exhaust  flow  expands  at  the  nozzle  lip  so  that 
the  pressures  of  exhaust  and  ambient  gases  are  matched  along  the  plume 
interface.  The  barrel  shock  forms  in  the  single  nozzle  plume  because 
expansion  waves  (upward  running  characteristics)  in  the  flow  reflect  from  the 
(near)  constant  pressure  plume  interface  resulting  in  reflected  compression 
waves.  These  eventually  focus  to  start  shock  system  ( B1 ) { see  Fig.  3-1 ( a ) ) . 

The  expansion  waves  which  start  this  process  can  arise  in  the  conical  like 
source  flow  leaving  the  nozzle,  however,  even  a  uniform  parallel  exit  flow 
nozzle  will  produce  the  same  result.  The  upward  running  characteristics 
leaving  the  exit  plane  become  expansion  waves  as  they  cross  the  Prandtl  Meyer 
expansion  fan  at  the  nozzle  lip  because  the  flow  is  axi symmetric  and  spreading 
laterally.  The  axi symmetric  nature  of  the  flow  causes  the  wave  strength  of 
the  B1  shock.  Figure  3-1 (b),  to  Increase  as  it  progresses  downstream  and 
approaches  the  axis  of  the  plume.  A  Mach  disc  and  reflected  shock  system  (B2) 
develop  downstream  of  the  original  barrel  shock  B1 .  The  flow  behind  the  Mach 
disc  is  subsonic  so  that  the  location  of  the  disc  depends  on  the  pressure 
distribution  and  mixing  processes  downstream  of  it.  Section  7  will  discuss 
this  portion  of  the  flow  field  in  detail.  This  is  in  distinction  to  the 
remainder  of  the  flow  which  is  supersonic  and  where  there  is  no  upstream 
influence.  This  inviscid  flow  pattern  is  well  understood  and  several  computer 
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Fig.  3-1  Single  Axisymmetric  Plume  Flow  Field 
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codes  are  available  (in  varying  degrees  of  approximation)  to  predict  it  (Ref. 
3-1,  3-2). 

It  is  quite  informative  to  investigate  the  nature  of  the  impingement  of 
two  uniform  jets  before  considering  the  complete  problem  of  the  impingement  of 
two  underexpanded  plumes.  This  flow  pattern,  outlined  schematically  in  Figure 
3-2  is  expected  to  have  two  shock  wave  systems.  The  impingement  shock  (I)  and 
the  recompression  shock  R.  In  the  side  view  the  I  shock  appears  basically  as 
expected  from  a  two  dimensional  pattern.  The  impingement  shock  turns  the 
oncoming  flow  parallel  to  the  symmetry  plane.  A  complex  process  takes  place 
at  the  intersection  of  the  plume  boundary  and  the  I  shock.  Based  on  work  by 
Hunt  and  coworkers  (Ref.  3-3,  3-4,  3-5)  the  discontinuous  boundary  pattern 
sketched  in  Fig.  3-2  is  expected.  These  references  deal  with  normal 
impingement  of  uniform  jets;  however,  the  interaction  of  the  I  shock  and  the 
plume  boundary  is  locally  equivalent  to  that  when  viewed  in  a  coordinate 
system  parallel  to  the  shock/boundary  intersection.  The  flow  field  is 
projected  onto  a  plane  perpendicular  to  the  intersection  line  (Fig.  3-2 ( b ) ) . 
The  component  of  velocity  parallel  to  the  intersection  line  is  constant  in  the 
neighborhood  of  the  intersection  line  because  (a)  it  is  parallel  to  the  shock 
wave  and  is  hence  unaltered  by  it  and  (b)  is  locally  parallel  to  the  plume 
boundary  both  upstream  and  downstream  of  the  impingement  shock.  The  pressure 
at  points  A  and  D  are  matched.  Downstream  of  the  impingement  shock  the 
pressure  at  B  is  greater  than  at  A  and  therefore  greater  than  at  D. 

Therefore,  an  expansion  fan  emanates  from  the  plume  boundary  at  the  point  of 
impingement  to  cancel  the  pressure  rise  due  to  the  I  shock  wave  (Station  1, 
Figure  3-2 (a)).  The  required  pressure  match  downstream  of  the  shock  demands 
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Fig.  3-2  Schematic  Flow  Field  for  the  Impingement  of  Two  Uniform  Plume*  (Sheet  1  of  3  ) 
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Fig.  3-2  Schematic  Flow  Field  for  the  Impingement  of  Two  Uniform  Plumes  (Sheet  2  of  3) 


that  at  the  point  of  impingement  there  is  at  least  sonic  velocity,  relative  to 
the  intersection  line,  to  support  a  Prandtl  Meyer  fan.  At  Station  2,  a  new 
feature  develops  in  the  flow-expansion  wave  fronts  stretching  in  three 
dimensions  interact  with  the  constant  pressure  boundary  giving  rise  to  inward 
moving  compression  wave  surfaces  that  coalesce  to  form  a  recompression  ( R ) 
shock  system.  A  schematic  of  this  detail  is  shown  in  Fig.  3-2  (c).  This 
coalescence  is  completely  analogous  to  the  formation  of  the  barrel  shock  (B) 
system  in  the  axi symmetric  case.  Subsequently,  (Station  3-5)  the  R  shock 
system  shrinks  in  size  and  grows  in  strength  as  it  approaches  the  plume 
center.  Another  way  of  viewing  the  overall  impingement  process  is  to  consider 
that  the  impingement  shock  by  elevating  the  pressure  of  a  perfectly  matched 
plume  creates  an  underexpanded  jet  which  subsequently  expands  laterally  giving 
rise  to  the  shock  pattern  familiar  to  underexpanded  plumes. 

The  flow  pattern  associated  with  the  impingement  of  two  underexpanded 
plumes,  in  general,  contains  all  three  shock  systems:  the  barrel  shock,  the 
impingement  shock  and  the  recompression  shock.  These  shock  surfaces  propagate 
across  the  plume  flow  fields,  interact  and  give  rise  to  subsequent  generations 
of  shock  surfaces.  It  becomes  harder  and  harder  to  identify  each  shock 
specifically  in  succeeding  generations.  Many  shock  configurations  are 
possible  depending  on  the  relative  strengths  of  the  three  systems  and  the 
order  in  which  they  intersect.  Three  observed  configurations  will  be 
discussed.  Each  flow  schematic  is  followed  by  a  corresponding  glow  photograph 
(Ref.  3-6). 

In  the  weakest  type  interaction  (Fig.  3-3)  the  flow  pattern  is  initially 
that  of  two  individual  plumes.  At  distances  less  than  the  first  Mach  discs 
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(C)  DETAIL  SHOWING  THE  FORMATION  OF  THE  RECOMPRESSION  SHOCK 
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the  plumes  appear  almost  as  two  individual  plumes  in  the  top  view.  The  first 
Mach  cell  is  only  slightly  distorted  by  the  I  shock  (see  top  view).  In  the 
top  view  the  impingement  shock  seems  to  merge  with  the  barrel  shock  and 
deflect  so  that  the  Mach  discs  are  not  quite  normal  to  the  nozzle  centerline. 
The  next  major  shock  pattern  occurs  downstream  of  the  Mach  discs  in  the 
central  portion  of  the  flow  between  the  exit  of  the  two  nozzles.  The  R  shocks 
(side  view.  Fig.  3-3)  from  the  upper  and  lower  portion  of  the  flow  intersect 
to  form  a  wedge  shaped  shock  pattern  in  the  flow.  The  leading  edge  of  this 
system  is  cut  off  (point  a  in  top  view  of  Figure  3-3)  as  it  is  intersected  by 
the  reflected  barrel  shock  downstream  of  the  Mach  disc.  The  small  features 
adjacent  to  this  central  pattern  (see  top  view  of  glow  photograph  Figure  3-3) 
seems  to  be  created  at  the  intersection  of  the  reflected  barrel  shock  and  the 
recompression  shock  at  point  b  shown  in  the  side  view  of  Figure  3-3. 

At  lower  background  pressures  the  initial  expansion  at  the  nozzle  lip  is 
greater  and  the  plumes  impinge  at  higher  angles  increasing  the  strength  of  the 
impingement  shock.  Figure  3-4  is  an  example  of  a  moderate  interaction  where 
the  impingement  shock  strength  is  increased  to  the  point  -where  it  cuts  off  the 
barrel  shock  system  before  the  formation  of  the  Mach  disc  associated  purely 
with  the  barrel  shock.  In  this  case  downstream  of  the  B/R  intersection  (top 
view)  the  R  and  transmitted  B  shock  intersect  in  such  a  way  as  to  create  a 
normal  shock  (Mach  disc)  in  the  center  of  the  flow.  In  the  strong  interaction 
case  (Fig.  3-5)  the  impingement  shock  rapidly  traverses  the  plume  and  diverts 
the  B  shock  sharply  toward  the  symmetry  plane.  This  transmitted  B  shock 
reaches  the  symmetry  plane  (top  view)  at  point  a  while  the  R  shock  (side  view) 
is  still  out  near  the  plume  boundary.  As  the  8  shock  system  reflects  from  the 


Fig.  3-4  Flow  Pattern  of  Underexpanded  Twin  Plumes,  Intermediate  Interaction 


symmetry  plane  a  V  shaped  trace  is  created  in  the  side  view  (Fig.  3-5). 
Subsequently  this  reflected  B  shock  intersects  with  the  R  shock  surfaces 
producing  an  irregularly  shaped  leading  edge  because  both  these  shock  surfaces 
are  not  planar  (Fig.  3-5,  station  2). 

The  three  basic  flow  patterns  described  most  probably  represent  only  a 
fraction  of  the  possible  flow  configurations.  There  are  likely  many 
variations  of  these  patterns  and  others  not  yet  observed.  For  this  reason  any 
computational  scheme  chosen  to  pursue  solutions  of  these  flow  fields  cannot  be 
of  the  type  that  is  constrained  by  geometrical  and/or  topological 
limitations.  The  scheme  must  allow  for  a  wide  variety  of  geometrical  patterns 
and  must  be  flexible  enough  to  include  as  yet  unknown  and  unanticipated 
configurations.  A  computational  technique  satisfying  this  constraint  is 
discussed  in  the  Section  5.  The  method  produces  accuracy  and  maintains 
efficiency  by  incorporating  detailed  local  flow  solutions  wherever  possible. 
Thus  shock  waves  and  slip  surfaces  are  portrayed  as  di scontinuties  which  are 
tracked  individually  and  locally  satisfy  the  appropriate  jump  conditions. 

There  are  other  three  dimensional  flow  features  in  these  plumes  that  must 
likewise  be  modelled  in  the  small  by  their  local  solutions.  Two  such 
solutions  required  for  the  multi  nozzle  plume  flow  field  have  not  been  analyzed 
in  the  past.  A  discussion  of  the  nature  of  the  lift  off  of  the  impingement 
shock  and  the  transition  of  the  regular  reflection  to  a  Mach  reflection 
process  in  the  intersection  of  two  three  dimensional  shock  surfaces  is 
presented  in  Section  6. 
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SECTION  4 


RELATIONSHIP  OF  IR  SIGNATURE  TO  PLUME  PROPERTIES 

There  is  a  variety  of  programs  which  require  knowledge  of  rocket  plume 
flow  fields.  Detailed  spatial  maps  of  exhaust  temperatures,  pressures  and 
chemical  species  concentrations  are  required  as  input  in  electromagnetic 
attenuation  codes,  aircraft  or  spacecraft  impingement  analysis  and  IR  emission 
codes  for  both  heat  transfer  and  optical  signature  evaluations.  These 
calculations  are  often  so  expensive  to  perform  that  parametric  analysis  or 
calculation  of  a  large  member  of  data  sets  is  out  of  reach.  For  that  reason 
it  is  difficult  to  answer  questions  pertaining  to  the  variation  of  plume  IR 
signature  as  a  result  of  systematically  changing  the  many  input  parameters. 
Normal  system  design  procedures  become  extremely  costly  or  have  to  be 
bypassed.  The  main  concern  in  this  section  of  the  study  was  the  effect  of 
multinozzle  plume  flow  fields  on  the  IR  signature  of  a  missile.  The  results 
were  presented  in  detail  in  Ref.  4-1  and  will  be  reviewed  briefly  here. 
Analytic  formulas  were  derived  which  relate  IR  signals  to  various  plume 
properties  thus  alleviating  to  some  extent  the  problems  in  parametric 
analysis.  These  relationships  point  out  the  allowable  errors  or  uncertainties 
in  plume  models  that  lead  to  desired  levels  of  accuracy  in  overall  IR 
predictions.  It  was  determined  that  there  are  several  other  processes  besides 
multi  nozzle  plume  impingement  that  must  be  properly  accounted  for  if  accurate 
IR  predictions  are  to  be  achieved  over  the  entire  altitude  range  .  The 
analysis  Is  best  suited  for  optically  thin  plumes,  however,  as  in  other 
situations  the  conclusions  will  probably  have  a  much  broader  range  of 
appl icabil i ty. 
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The  underlying  source  term  in  all  optical  radiation  calculations  is  the 


Planck  Black  Body  Function.  The  radiant  emission  from  elemental  gaseous 
sources  is  proportional  to  the  product  of  that  function  and  the  absorption 
coefficient  of  the  IR  active  molecule.  The  magnitude  of  the  Planck  function 
is  sensitive  to  temperature  and  becomes  increasingly  so  as  the  temperature 
decreases.  This  sensitivity  is  qreater  the  shorter  the  wavelength  as  shown  in 
Fig.  4-1  where  the  logarithm  (base  10)  of  the  Planck  function  multiplied  byA.^ 
is  plotted  versus  temperature.  Temperature  sensitivity  is  graphically 
depicted  in  Fig.  4-2  where  the  temperature  increment  necessary  to  produce  an 
increase  in  the  Planck  function  by  a  ratio  of  2.71,  1.65,  and  1.28  (e,  e^, 
e^4)  is  shown  for  the  range  400-2000  K.  In  Fig.  4-2a,  for  example,  for  2.7 
microns  at  1000  K,  a  temperature  increment  of  100  K  produces  an  increase  in 
the  Planck  function  by  a  factor  1.65;  at  2000  K,  the  same  increase  requires  a 
temperature  increment  of  380  K.  In  fact,  the  same  accuracy  at  800  K  requires 
a  temperature  accuracy  of  60  K.  Figures  4-2b  and  4-2c  show  the  temperature 
sensitivity  for  shorter  and  longer  wavelengths,  and,  as  before,  the  decreasing 
sensitivity  with  increasing  wavelength  is  shown.  Below  the  abscissa  in  Fig. 
4-2a  there  is  an  approximate  altitude  corresponding  to  the  temperature  axis. 
This  correspondence  is  approximate  but  serves  as  a  reasonable  guide.  It  is 
interesting  to  note  that  uncertainties,  errors  in  calculation  or  variations  in 
any  plume  parameter  that  lead  to  temperature  increments  have  much  less  of  an 
effect  at  low  altitudes.  Thus  a  change  in  temperature  at  sea  level  of  400  K 
(all  other  quantities  fixed)  would  only  produce  a  change  of  the  Planck 
function  by  a  factor  of  1.65.  On  the  other  hand  at  50  km  altitude  there  is  a 
factor  of  2.71  for  only  150  K  variation.  This  critical  temperature  dependence 
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Fig.  4-1  Planck  Black  Body  Function  Variation  with  Temperature 
for  Several  I R  Wave-lengths 


requires  plume  models  to  have  very  tiqht  controls  on  the  approximations 
involved  in  the  prediction  of  plume  temperature  levels  in  order  to  achieve 
accuracy  over  the  entire  altitude  range. 

A  variety  of  processes  control  the  temperature  levels  in  the  plume 
gases;  expansion  to  background  pressure,  turbulent  mixing  of  the  exhaust  and 
ambient  gases,  and  afterburning  cnemistry  in  the  mixing  between  the  exhaust 
gas  and  ambient  stream.  The  mixing  process  is  of  direct  concern  because  it: 
adds  a  temperature  increment  due  to  viscous  dissipation  and  controls  the 
geometric  size  of  the  radiating  region.  At  low  altitudes,  an  equally 
important  contribution  to  the  radiation  levels  is  the  afterburning  chemistry, 
which  not  only  adds  a  temperature  increment,  but  can  also  be  responsible  for 
the  consumption  or  production  of  a  radiating  species.  The  influence  of  both 
mixing  and  afterburning  has  been  the  target  of  other  studies  and  will  not  be 
focussed  on  here.  The  inviscid  expansion  process  of  the  exhaust  gas  to  the 
ambient  pressure  levels  becomes  more  and  more  significant  as  altitude 
increases.  There  are  two  principle  processes  that  cause  major  changes  to  the 
temperature  that  would  oe  achieved  via  an  isentropic  inviscid  expansion  from 
the  nozzle  exit  to  the  background  pressure.  The  primary  deviation  is  due  to 
shock  waves  in  the  plume  flow  which  cause  large  entropy  increases  that  Ders’st 
into  the  plume  far  field.  These  shock  waves  result  from  the  adjustment  of  the 
underexpanded  nozzle  flow  to  the  local  pressure  of  the  surrounding  fluid  ana 
the  impingement  of  individual  exhaust  plumes  on  each  other  in  the  case  where 
the  vehicle  has  more  than  a  single  engine.  The  former  case  is  well  understood 
because  it  is  an  axisymmetric  flow,  while  the  latter  is  a  complex 
three-dimensional  flow  that  was  the  object  of  this  research  Drogran. 
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In  order  to  properly  account  for  large  pressure  ratio  expansions  in 
multinozzle  plumes  all  flow  processes  beyond  shock  heating  that  influence  the 
temoerature  field  should  be  examined.  The  role  of  nonideal  thermodynamic 
properties  of  the  fluid  in  the  inviscid  expansion  has  to  be  properly  addressed 
in  any  plume  model.  The  temperature  difference  between  ideal  and  nonideal 
expansion  increases  with  overall  pressure  ratio  and  hence  altitude.  The 
nonideal  expansion  process  reduces  the  plume  core  temperature  and  therefore 
further  compounds  temperature  sensitivity.  Two  common  assumptions  that  must 
be  reviewed  are:  (1)  the  exhaust  medium  is  calorically  perfect  and  (2)  the 
exhaust  composition  is  frozen.  Exhaust  gases  can  contain  substantial  mole 
fractions  of  water  vapor  and/or  carbon  dioxide  that  exhibit  changing  values  of 
specific  heat  over  the  entire  temperature  ranqe  (500-2000  K).  Other  triatomic 
molecules  can  exhibit  similar  specific  heat  temperature  dependence  over  these 
temperature  ranges.  This  results  in  a  substantial  temperature  decrease  due  to 
the  difference  between  a  frozen  composition  and  the  more  widely  employed 
constant  V  expansion.  In  addition  all  finite  rate  chemical  reactions  cannot  oe 
considered  frozen.  Thus  changing  composition  must  be  examined. 

Assessment  of  the  gross  effects  of  multiple  nozzle  self  impinqement  ann 
the  real  gas  expansion  process  on  the  plume  IR  signature  was  achieved  v’a  an 
analytic  analysis  derived  and  described  in  detail  *n  Ref.  4-1.  The  analyses 
is  for  optically  thin  conditions  and  as  such  does  not  account  to r  se'f 
absorption  of  the  outer  cooler  regions  of  the  olume.  'he  mathematical 
analysis  for  the  far  field  is  based  on  standard  methods  which  take  advantaae 
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of  the  exponential  nature  of  the  Planck  function.  The  results  of  the  analysis 
was  tested  against  computer  generated  plume  signatures.  Two  model  plumes  were 
generated  and  the  resulting  station  radiation  was  computed  using  GRUMPLUME 
(Ref  4-2).  These  results  were  used  to  show  that  the  computed  station 
radiation  was  proportional  to  the  derived  formulas. 

A  schematic  of  the  model  plume  flow  field  is  shown  in  Fig.  4-3.  The  near 
field  multinozzle  impingement  and  the  subsequent  expansion  to  amoient  pressure 
are  assumed  to  have  taken  place  upstream.  The  model  studied  here  is  a  far 
field  model  and  does  not  address  the  emission  from  the  regions  near  the  exit 
plane.  The  initial  impingement  region  is  highly  nonuniform  and  requires  a 
much  more  complex  formulation  that  would  not  be  expected  to  yield  simple 
relationships.  Two  plumes  were  investigated  corresponding  to  altitudes  of  50 
and  60  km  with  start  line  temperatures  of  383  K  and  745  K  respectively  and 
water  vapor  mole  fractions  of  .33  and  .32.  The  computed  centerline 
properties  from  GRUMPLUME  are  shown  in  Fig.  4-4. 

The  station  radiation  predicted  for  these  two  test  cases  is  shown 
plotted  in  a  normalized  form  in  Fig.  4-5.  The  striking  feature  of  the  result 
is  that  in  the  initial  region  of  the  far  field  the  station  radiation  is  a 
linear  function  of  distance.  The  results  of  the  analysis  are  as  follows.  In 
the  initial  region  of  the  far  fi el d  the  initial  value  of  the  station  radiat’on 
at  z  =  0  is  proportional  to  the  product  of  the  2lanck  function  at  the 
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Fig.  4-3  Schematic  of  Model  Plume  Flow  Field 
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Fig.  4-4  Centerline  Temperature  and  Water  Vapor  Mole  Fraction 
for  Model  Plume* 
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Fig.  4-5  Normalized  Station  Radiation  for  the  Two  Model  Plumes 


Fig.  4-6  Correlation  of  Far  Field  Station  Radiation 
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initial  temperature  and  the  initial  value  of  the  mole  fraction  of  water  vapor 
(the  radiating  species)  .  This  formula  correlated  the  results  very  well  as 
shown  in  the  table  below. 


h 

- 1 - 

J2(0) 

Jz(0) 

1 

1 

i 

t 

XH2o(°)  nx  °(°) 

...  _ 

50 

2.678 

3390 

60 

|  .828 

3310 

(a! 1  units  arbitrary ) 

The  result  is  only  strictly  applicable  for  the  distribution  given  by  the  model 
plume.  In  a  more  accurate  fluid  dynamic  model  of  the  plume  the  initial 
temperature  profile  would  not  be  constant  resulting  in  a  more  complex 
relationship.  In  the  plume  far  field,  where  the  centerline  temperature  and 
species  mole  fraction  of  water  vapor  decay  both  axially  and  radially  the 
results  are  (see  Ref.  4-1  for  details  of  the  derivation) 

Jz  a  NA°(Tq)  pT,- /(d2I/dr2)^ 

The  station  radiation  is  proportional  to  the  product  of  mole  fraction  of 
radiating  species,  Planck  function,  static  pressure  and  temperature  and 
divided  by  the  second  derivative  of  temperature  all  evaluated  at  centerline 
conditions.  The  validity  of  this  correlation  is  demonstrated  in  Fig.  4-6 
where  the  entire  series  of  points  in  the  far  field  of  both  plumes  are  plotted 
versus  centerline  temperature.  All  the  points  fall  within  approximately  5»  of 
the  value  1 .45. 
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Two  plumes  are  sketched  in  Fig.  4-7:  a  single  engine  axisymmetric  plume 


and  a  multiple  engine  three-dimensional  plume.  The  shock  heatinq  in  the 
axisymmetric  case  gives  rise  to  a  very  highly  peaked  temperature  profile 
because  the  initial  shock  strength  is  weak  and  increases  rapidly  as  it  nears 
the  axis.  The  fraction  of  plume  mass  flow  the  shock  intercepts  when  it  is 
strong  is  quite  low  because  the  streamlines  are  continuously  diverging  from 
the  axis  in  the  plume  core.  In  the  multiple  engine  case,  a  strong  shock 
caused  by  the  plume  impingement  intercepts  a  large  fraction  of  the  exhaust 
flow  near  the  exit  plane.  Therefore,  the  temperature  profile  downstream  will 
have  a  flat  shape  and  a  larger  fraction  of  the  exhaust  mass  flow  will  be 
heated  to  high  temperatures.  In  the  analysis  presented  the  effect  of 
oercentage  mass  flow  affected  by  shock  heating  was  not  explicitly  addressed, 
however,  the  station  radiation  will  be  proportional  to  that  percentage.  In 
the  previous  discussions  it  has  been  established  that  the  station  radiation  is 
directly  proportional  to  the  Planck  function  based  on  the  peak  temperature. 
Figure  4-8  demonstrates  what  can  be  expected  to  be  the  variation  of  Planck 
function  due  to  shock  heating.  The  ordinate  is  the  ratio  of  Planck  function 
at  shock  heated  temperature  to  Planck  function  at  the  isentropic  expansion 
temperature.  The  abscissa  is  the  isentropic  temperature  achieved  by  an 
expansion  from  2000  K.  Below  this  is  an  additional  scale  showing -the 
corresponding  altitude.  Figure  4-3  includes  a  series  of  curves  for  an  entire 
range  of  total  pressure  losses.  The  greater  the  total  pressure  loss  the 
higher  the  temperature  and  hence  the  higher  the  ordinate.  Notice  the  size  of 
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Fig.  4-7  Schematic  of  Single  Engine  &  Multiengine  Exhaust  Plumas 
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a  factor  of  2  shown  on  the  chart.  It  is  clear  from  this  chart  that  multiple 
nozzle  impingement  effects  become  progressively  more  influential  as  altitude 
increases.  At  low  altitidues,  below  say  10-20  km,  the  individual  plumes 
interact  only  weakly  so  that  there  is  little  additional  shock  heating 
developed  and  Fig.  4-8  shows  that  shock  heating  has  diminishing  influence  at 
those  altitudes  in  any  case. 

The  sensitivity  of  the  Planck  function  to  temperature  variations 

requires  tight  controls  on  all  processes  both  numerical  and  physical  that  lead 

to  error  beyond  the  shock  heating  caused  by  multiple  nozzle  self  impingement. 

There  are  a  variety  of  considerations  that  are  discussed  in  some  detail  in 

Ref.  4-1.  To  complete  the  discussion  in  this  section  the  role  of  nonideal 

thermochemi cal  properties  will  be  examined.  In  Fig.  4-9  the  temperature 

achieved  by  expansion  of  a  typical  exhaust  gas  composition  to  ambient  pressure 

is  plotted.  Three  isentropic  expansions  are  shown.  Two  are  for  ideal  gases 

with?=  1.22  and  7  =  1.3  (y  is  the  isentropic  exponent  C  /C  )  and  one  for 

P  v 

frozen  composition.  The  nozzle  exit  plane  value  of  Y  “is  1.22.  It  is  evident 
from  Fig.  4-9  that  above  an  altitude  of  about  20  km  the  assumption  of  a 
constant  value  of  y  will  overpredict  exhaust  gas  temperatures.  The  frozen 
composition  expansion  passes  smoothly  between  the  bounds  provided  by  the  two 
constant  7  expansions.  The  frozen  chemical  composition  expansion  process  is, 
therefore,  a  more  desirable  feature  to  employ  in  an  accurate  plume  model.  The 
accuracy  of  the  frozen  expansion  process  must  be  examined  in  the  light  of  the 
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Fig.  4*9  Comparison  of  Frozan  Composition  and  Constant  7  Expansions 


fact  that  chemical  reactions  do  not  necessarily  freeze  out  even  at  higher 
altitudes.  Figure  4-10  (Ref.  4-1)  displays  the  Damkohler  number  as  a  function 
of  altitude  for  the  reaction  set  shown.  Reactions  4,  5,  6,  9  can  still  be 
active  at  the  high  end  of  the  altitude  range.  The  possible  shuffling  of 
chemical  species  by  these  four  reactions  must  be  examined  to  determine  if  they 
can  substantially  alter  the  temperatures  achieved  by  the  frozen  expansion. 
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Tab!#  I  Typical  Reaction  Set  For  Amine  Fuel  Exhaust 
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Fig.  4*10  Damkohler  Number  for  RMction 
Sot  ai  a  Function  of  Altitude 
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SECTION  5 


NUMERICAL  CALCULATION  PROCEDURE 

The  variety  of  possible  structures  of  the  multiple  nozzle  plume 
impingement  flow  field  requi res  that  the  computational  procedure  employed  in 
their  analysis  be  as  independent  of  geometry  as  possible.  The  code  must  allow 
for  a  wide  variety  of  shock  configurations  and  have  the  flexibility  to  allow 
for  as  yet  unknown  additional  geometries.  For  this  reason  the  "floating 
discrete  shock  fitting"  approach  originally  devised  by  Moretti  was  chosen 
(Reference  5-1,  5-2).  A  two  dimensional  version  of  this  method  was  employed 
successfully  by  Salas  (Reference  5-3)  for  a  complex  two  dimensional  scramjet 
flow  field  containing  many  shock  waves.  For  the  present  problem  these  methods 
were  generalized  to  compute  three  dimensional  steady  inviscid  flow.  Three 
types  of  mesh  points  are  recognized:  interior,  boundary  and  discontinuity,  The 
computational  mesh  is  a  fixed  Cartesian  grid  where  the  shock  surfaces  and 
pressure  boundaries  propagate  freely  across  the  grid.  The  shock  and  pressure 
boundary  surfaces  are  portrayed  as  discontinuity  surfaces  and,  for  example, 
the  complex  impingement  shock/boundary  interaction  is  modelled  in  detail  as  a 
point  singularity. 

The  discontinuity  surfaces  are  considered  on  a  cell  by  cell  basis  so  that 
the  computer  code  must  contain  the  necessary  logic  to  oerform  the  correct 
calculations  in  all  possible  cases.  Because  the  cross  angle  of  the 
discontinuity  surfaces  must  be  determined  to  compute  propagation  velocity  the 
discontinuity  cells  must  be  connected.  Thus  the  computations  necessary  for 
each  discontinuity  cell  cannot  be  done  independent  of  properties  in 
neighboring  cells.  The  number  of  possible  configurations  is  large  but 
manageable  for  single  snock  surfaces  because  it  is  possible  to  devise  rules  to 
cover  a  large  portion  of  the  possibilities  and  to  identify  and  code  for  the 


exceptions.  When  two  shock  surfaces  intersect. 


tne  geometric  oossibi ’ i ties 


and  hence  complexity  increases  enonmously  witn  the  attendant  'ncreases  ;n 
program  logic.  The  marching  step  size  is  restricted  to  maintain  the  same 
configuration  of  discontinuity  cells  for  each  step.  At  the  end  of  eacn  step, 
if  any  of  the  discontinuities  has  reached  a  mesh  point,  it  is  crossed  to  the 
other  side  and  all  the  necessary  indicators  wnich  guide  the  program  logic  are 
reset  to  account  for  the  new  configuration. 

The  exhaust  gas  is  considered  inviscid,  thermally  and  calorically  perfect 
and  is  governed  by  the  three  dimensional  Euler  equations.  In  Cartesian 
coordinates,  the  conservation  of  mass,  momentum  and  entropy  are 


ux  *  vy  w2  +  y  (uPx  +  vPy  +  wPz)  *  0  (5-1.1) 

UUX  +  vUy  +  WUZ  +  TPX  =0  (5-1.?) 

uvx  +■  Wy  +  wvz  +  TPy  =  0  (5-1.3) 

uwx  +  vWy  +  wwz  +  TPZ  -  0  (5-1. A) 

uSx  +  vSy  +  wSz  =  0  (5-1.5) 


where  u,  v,  w  are  the  Cartesian  velocity  components  in  the  x,  y ,  z  directions 
and  P  is  the  natural  logarithm  of  the  pressure,  S  the  entropy,  T  the 
temperature,  7  the  ratio  of  specific  heats  and  all  thermodynamic  quantities 
are  non  di mens ional i zed  by  the  stagnation  conditions  (e.g.,  o0,  P0,T0), 
the  velocity  components  with  respect  to  the  quanti ty^^/sypT  and  the 
coordinates  by  the  nozzle  exit  radius.  The  marching  (time  like)  direction  is 
the  z  direction  and  the  following  two  flow  angles  are  introduced 


3'C 


f  =  u,  w  5-  1 .  i 

(7=  7  X  5  -  1  .  ' 

Suotract  £a.  (5-1. 4  mult  'Ql '  ed  by  j  *,'om  -  j,  5-1.-  nu't*2'-ed  ov  v  jrc 
divide  the  results  by  w-  to  get 

m 

rr  ♦  at  -  r  -  7  w“  ?  -  rP  )  3  3  • 

x  j  z  <  2  '  -  *  * 

Similarly  an  equation  *or  a  i s,  der ’ ved  *-om  Eqs.  5-1.4  and  5-1.5 

rax  tfO’y  r  J,  7  w-;  3/  -  s3,.  =  5  5-5.5 

Subtract  Eq.  (5-1.4)  from  Eo.  '5-1.1'  t ' mes  h  ana  qivde  tne  -esu't  ov  *2 

to  get 

■m 

Tx  +  a^  +  V(TPX  *<7P  +  (l-7,yw)  '■  =  5  '5-2.2' 

Divide  Eq.  (5-1.5)  by  w  to  get 

rS  +  <?S  +  S  =  0  5-?.  4 

x  J  z 

Equations  (5-2.1  -  5-2.4)  form  a  set  of  four  equations  for  tne  four  o r-marv 
unknowns  r  ,  a  ,  3  and  S.  The  remainder  of  tne  ow  quantities  are  :er-  /ea 
from  the  following  subsidiary  equations.  Employing  tne  donservat 1on  if 
stagnation  enthaloy  the  equation  for  w  (the  axial  ve'ocity 

?  2Y(l-7) 

■>-l)  (l+t^-HT)  1-  1 

The  equation  of  state  for  7 

7  =  exp  (((y-l)P  +  s)/y)  5-2.5 

7he  definitions  (5-1.6)  and  (5-1.7)  for  u,  v 

u  -  rw  5-2.3' 

v  a  aw 
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’  ne  nter’ jr  301  its  ere  calculated  dv  tie  integration  of  Eos.  (5-2.1  -  5-2. 4) 
emo'ov.nq  tie  two  steo  ‘‘'acCormacx  scneme  (3er‘.  5-4).  "he  marcninq  direction  is 
tie  z  axis  and  tne  sequence  of  ca" tu'ations  1  s  to  solve  **or  ^rom  Ed. 

5-2.3',,  tnen  vs’ig  tms  value  -n  (5-2.1:  and  (5-2.2)  solve  for  rz  ana  <j7. 

3Z  is  comouted  di^ectlv  ‘>om  Ed.  '5-2. 4).  "ne  oasic  notation  for  mesh 
oomts  and  :e:'s  's  snown  •  n  :’j.  5-1,  .vnicn  oeoicts  a  segment  of  the 
tomoutat ’onal  o'ane.  Eac.n  ”esn  omit  s  'ocateo  ov  a  single  integer  ana 
dent  1  *  ’  ed  ,r>  *ne  'igure  ov  an  jooer  case  'etter.  A'ong  each  row  the  mesn 
joint  iumoer  iic-eases  one  oy  nov’iq  one  mesn  001  it  to  the  right  and  increases 
oy  3,  tne  number  of  ooi  nts  'i  eacn  '■ow,  oy  movincj  jd  one  mesn  ooint.  Each 
ce1 1  -s  ’lentif’ed  oy  tie  /a ‘ oe  of  J  of  tne  lower  left  hand  mesn  ooint  and 
denoted  oy  "‘ms  s’nqie  an  ay  notation,  though  unconventional  ,  orocuces 
some  orogramm-ig  s 1 mp ; • * ’cat’cn ,  ano  does  not  increase  computer  operations 
s’rce  ne  norma'  icub  e  sunsniDt  approach  reauires  the  transparent  to  the 
,ser  same  computafons  to  '  ocat  ?  3  variable. 

*he  como  ete  ce,’a’’s  of  the  use  of  the  ^acCormacK  finite  aif*erence 
3  qor’tnm  re  y  •  n  many  o’aces  3ef .  5-4)  and  w’ ' '  not  oe  -eoeatad 

mere.  In  tie  'itegrat'on  scheme,  f3nsverse  oart’al  terivat’/es  'der- vat' ves 
1  “ne  o'ane  normal  to  tne  marr.i’rg  01 -ect  1  on  are  evaluated  Oy  two  ooint. 

*' 0  rwa  rq  or  oacxwarq  ‘jrmu’as.  "he  not  at' on  for  -ne  transverse  derivatives  in 
snown  -n  5  <  a.  5-2.  -nr  any  y  tne  **'  ow  tua  nt  ■  ’ '  as  ,  tne  fo  1 '  owi  na  'omu’a  are 
'■'n”.e  ”,‘3r0f'C3S  at  oomt  J. 
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Fig.  5-2  Notation  for  Forward  and  Backward  Differences 
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Fig,  5-3  Sample  Finite  Difference  Mesh 
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where  V  is  any  flow  quantity  and  the  aquations  define  the  symbols  DFX,  DFV, 

DBX  and  D8Y.  Combinations  of  forword  and  backward  differences  for  the  x  and  y 
directions  is  permitted  in  the  predictor  step  so  long  as  it  is  -eversed  in 
corrector  step.  If  backward  differences  are  employed  for  both  x  and  y  in  the 
oredictor  step  and  forward  differences  in  the  corrector  step  the  new  /a lues  at 
Z  +  A  Z  at  point  J  depend  only  on  tne  seven  mesh  points  bounding  the  two 
shaded  cells  in  Fig.  5-2.  If  there  is  a  discontinuity,  a  shock  wave  for 
example,  in  either  of  these  shaded  cells,  then  the  computation  at  noint  J  must 
be  modified  to  account  for  it.  The  orientation  of  this  'finite  difference 
molecule"  dictates  a  computational  logic  in  terms  of  defining  how  mesh  points 
are  affected  by  discontinuity  cells. 

Figure  5-3  shows  a  typical  computational  mesh  denoting  the  mesh 
numbering  system  and  snowing  a  single  continuous  shocx  or  discontinuity 
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surface.  Three  types  of  points  are  indicated  on  the  figure:  interior  points, 
oounaary  points  and  discontinuity  ooints.  Boundary  points  3re  those  interior 
ooints  associated  with  a  finite  difference  molecule  'see  rig.  5-2;  whfcn  nave 
a  discontinuity  passing  througn  it  and  so  cannot  Pe  comoutec  employing  tne 
standard  interior  point  calculation  scheme.  The  steps  involved  in  one 
computational  step  are  outlined  in  Fig.  5-4.  The  necessity  *'or  step  ?  is  to 
maintain  the  overall  configuration,  shat  is,  ail  interior  points  remain 
interior  points,  boundary  points  remain  boundary  points  ,:or  noth  the  predictor 
and  corrector  portions  of  the  calculation. 

In  step  4  (Figure  5-4),  the  computation  of  the  properties  at  the 
discontinuity  and  boundary  points  of  each  cell  containing  a  discontinuity  is 
considered  individually.  The  discontinuity  surfaces  are  considered  to  oe 
oriented  so  as  to  recognize  the  high  and  low  pressure  sides.  The  projection 
of  the  discontinuity  surfaces  on  the  transverse  or  computational  plane  results 
in  a  continuous  line  (c.f..  Figure  5-3).  If  this  line  is  followed  from  a 
prescribed  starting  point,  for  aacn  cell,  the  discontinuity  can  enter  on  any 
one  of  four  faces  and  leave  on  any  of  the  remaining  three  faces  -  giving  a 
total  of  twelve  oossible  configurations  shown  in  Fig.  3-5,  In  the  computer 
program,  for  each  discontinuity  cell,  tne  properties  at  the  discontinuity 
point  which  is  the  entry  ooint  to  the  cell  are  associated  with  that  cell.  co r 
each  of  tne  twelve  possible  configurations  for  the  di scant 1 nu' ty  ceils  tne 
computer  code  contains  the  necessary  loqic  to  determine  which  mesh  ooints  ar? 
boundary  points  and  which  discontinuity  Points  are  to  be  ca'culated.  it  any 
step  the  code  contains  a  tab’s  of  tne  di  sr.ont  i  mji  tv  :e  s  a  no  the  program 
logic  determines  tne  method  of  :a  cu'at’on  -pr  -ao  or  *he  .n*ocwn  points.  A 
typical  example  of  a  di  scont’  nui  ty  ce  ;  m,.«n  r  5-5.  a  ;ooc< 
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Fig.  5-4  Flow  Diagram  for  Ona  Step  Using  Floating  Discontinuity 
Finita  Difference  Computer  Program 
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DISCONTINUITY 


enters  from  the  bottom  face  of  the  call  and  leaves  through  the  too  face.  To 
the  right  of  the  shock  the  flow  is  known  (i.e.  undisturbed  free  stream  or 
undisturbed  olume).  The  code  logic  determines,  for  this  configuration,  t.oat 
the  two  points  to  be  calculated  associated  with  the  j  cell  are  mesh  point  J 
and  shock  point  Sj.  The  unknowns  at  Sj  are  the  values  of  r,  a,  P  ana  S 
on  the  high  pressure  side  of  the  shock  and  the  shock  point  velocity  dfj/dz 
where  fj  is  the  fractional  distance  of  Sj  from  J  (see  Fig.  5-6). 

The  shock  configuration  shown  in  Figure  5-6  is  the  simplest 
possible.  Each  of  the  discontinuity  cells  is  a  type  1  as  defined  in  Figure 
5-5.  The  computation  of  the  shock  cross  flow  angle  is  computed  by  locating 
the  shock  points  on  either  side  of  the  point  of  interest.  This  is  achieved  by 
looking  ahead  and  behind  the  discontinuity  cell  j  and  locating  the  snoci< 
points  in  cell  j-R  and  j+R.  In  this  case,  the  location  of  the  shock  point  in 
j+R  (j-R)  is  at  y  =  Ay  (y  =  -Ay)  and  x  =  fj.pA*  (x  =  fj.pAx)  (x,  v 
both  measured  relative  to  the  mesh  ooint  J).  In  other  conf i gurations  where 
cell  j-R  is  type  4  or  type  11  (see  figure  5-5)  the  computation  of  the 
location  of  the  preceeding  snock  ooint  follows  a  different  formula.  Lixewise 
for  the  following  point.  In  the  configuration  shown,  the  ooint  J  is  a 
boundary  point  and  must  be  calculated  accordingly,.  For  each  triplet  of  types 
of  discontinuity  points  (in  Fiaure  5-6  this  is  a  1-1-1)  a  logic  pattern  is 
designed  to  compute  the  shock  cross  angle  correctly  and  to  determine  wnich 
point  is  a  boundary  ooint. 

"he  calculation  procedure  used  for  discontinuity  points  and 
boundary  points  is  outlined  below  for  a  two  dimensional  flow-the  extension  of 
these  methods  to  tnree  dimensions  is  straight  forward  geometrically  but 
reauires  some  analysis  pertinent  to  tnree  dimensional  character! sti c  theorv 
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and  will  be  discussed  in  subsequent  oaragraphs.  Three  ideas  must  oe  developed 
in  order  to  explain  this  procedure:  (1)  computation  of  characteristic 
relationships  using  finite  difference  algorithms  (2)  slant  marching  steps  and 
(3)  pseudo  point  values.  Points  A,  S  and  C  form  the  finite  difference 
molecule  for  the  computation  of  the  values  at  point  N  in  a  two  dimensional 
MacCormack  finite  difference  scheme  (Fig.  5-7).  (Note  that  the  intermediate 
values  employed  in  the  predictor/corrector  scheme  are  only  rotational 
simplifications  and  that  only  values  of  properties  at  A,  B  and  C  are  employed 
in  the  computation  of  N.)  The  object  here  is  to  show  how  the  values  computed 
by  the  finite  difference  calculation  can  be  used  to  compute  the  characteristic 
relationships  at  N  without  the  need  for  interpolati ng  the  values  at  the 
points  denoted  +  and  -  which  are  at  the  foot  of  the  characteristics  (Figure 
5-7).  Denote  the  velues  at  N  computed  by  a  second  order  finite  difference 
scheme  to  be  upg,  wcq,  Pprg,  Spr)  from  which  ‘'jtq  (  *  is  the  Prandtl 
Meyer  angle)  can  be  computed.  The  two  characteristic  rel ationshiDS  associated 
with  point  N  (second  order  accurate)  are 


downward  wave 

&8 +&v=  r 

ii 

kA 

1 

> 

upward  wave 

A0- 

v  fr  =  '*(  r_  *  rN) 
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MARCHING  DIRECTION 


Fig.  5-7  Finite  Difference  Molecule  for  Two  Dimensional  Flow 


Where  T  is  the  inhomogeneous  or  forcing  term  containing  velocity 
gradients  out  of  the  plane  (c.f.  Fef.  5-5  for  derivation  of  these 
equations).  Employing  the  notation  of  Fig.  5-7,  these  relationships  Decome 

+  =T>Z  +9-  + 

^ 

In  the  strai gntforward  aoplication  of  the  method  of  characteristics,  these  two 
equations  are  employed  to  solve  for  the  two  unknowns  9y  and  „  .  .  This 
requires  the  interpolation  of  data  to  determine  and  r  .  3ecause 

the  finite  difference  result  ana  the  character!  st -c  -esu't  must  oe  the  same  to 
second  order 


=  T  AZ 


tan  * 


=  f\Az  -  9^ 


Thus,  the  simDle  result  at  point  W  'or  the  two  characteri stic  relationshio  is 


*M  *  N  *  tan_i(drn/w,n;  +  Vd 

•  •j  r  j 


^  =  tan"1 (up^/w-^) -  v 


FD 


which  does  not  naqyjre  any  interpolation  or  any  scheme  other  than  the 


ordinary  finite  difference  algorithm  used  at  any  interior  Doint. 


Figure  5-8  shows  that  (for  the  discontinuity  point)  the  new  values 
desired  at  z  +  Az  do  not  lie  on  the  same  value  of  x  as  point  B,  rather  the 
point  N  is  elevated  by  \Ax.  The  computation  of  this  point  is  made  simply  by 
forming  a  directional  derivative  from  the  governing  equations.  For  any 
system  of  partial  differential  equations 


au  =  A  au 

az  ax 


+  B 


where  A  and  B  are  the  matrices  appropriate  to  the  equations  of  interest  and  U 
is  the  unknown  vector.  To  compute  the  values  at  point  N  use  the  system 


dU 

dz 


au  +  au  Ax 

az  axAz 


+  1 


Ax\au 

12  jax 


+  B 


where  I  is  the  identity  matrix  and  dU/dz  is  now  the  derivative  used  to  compute 
point  N  from  point  B  in  Fig.  5-8  by  a  "slanted  step."  The  finite  difference 
scheme  is  applied  to  the  right  hand  side  of  this  equation  using  exactly  the 
same  rules  as  are  used  at  any  interior  point  resulting  in  a  second  order 
accurate  predication  at  N. 

The  final  ingredient  necessary  for  the  boundary  point  calculation  is 
the  pseudo  or  projected  point.  Figure  5-9a  shows  the  typical  configuration 
encountered  in  the  "shock  between  the  mesh  points"  type  calculation.  In 
general  at  station  z  the  required  properties  are  known  at  mesh  points  A  and  B 
and  on  the  high  pressure  side  of  S  the  shock  point.  These  points  do  not  form 
the  standard  equally  spaced  finite  difference  molecule.  In  order  to  use  the 
same  algorithm  that  is  used  at  all  other  mesh  points,  the  finite  difference 
molecule  shown  in  Fig.  5- 9b  is  constructed.  The  values  of  the  flow  variables 
at  the  point  P  are  calculated  by  a  simple  linear  extrapolation  from  the  values 
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Fig.  5-9  Pseudo  Point  Construction 
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at  A  and  S.  The  value  at  3  is  not  used,  This  technique  nas  the  critical 

* 

oroperty  that  for  geometries  whp'-e  point  S  is  very  close  to  Doint  3  truncation 
errors  are  not  amplified  in  creating  the  value  at  This  molecule  can  oe 
shown  to  produce  a  first  order  rendition  of  the  distribution  of  flow 
properties  on  the  mesh  segment  ABS.  However,  a  simple  comparison  between  the 
Taylor  series  expansions  employing  properties  at  A,  B  and  S  with  one  employing 
A,  B  and  P  shows  that  the  error  involved  is,  in  the  worse  case, 

(1/8)  (d2f/  dx2)h2  (f  being  any  flow  property).  Thus  '•easonably  small 
errors  are  expected  from  this  approximation. 

The  sequence  of  steps  for  calculating  a  typical  boundary  and  shock  point 
(the  flow  upstream  of  the  shock  is  known  in  this  case)  are  as  follows  (refer 
to  Fig.  5-10):  (1)  the  shock  location  at  z+  A z  (SN)  is  estimated  based  on  the 
known  shock  slope  at  S;  (2)  the  flow  variables  at  point  P  are  computed:  (3) 
the  finite  difference  molecule  ABP  is  employed  to  compute  the  values  at  point 
N  using  the  standard  MacCormack  scheme;  (4)  the  slant,  step  and  same  procedure 
is  used  to  compute  the  finite  difference  values  of  the  new  variables  at  SN ; ( 5 ) 
these  new  values  are  combined  to  yield  a  single  characteristic  eauation  -  in 
this  case 

%N  "  "SN  =  tan  1  'uFD/wFD^“  '  FD  ’ 

(6)  the  Ranki ne-Hugoni ot  equations  and  two  estimates  of  the  shock  slope  at  S N 
are  used  to  calculate  two  solutions  for  the  conditions  behind  the  shock  at  SN: 

(7)  assume  linear  variation  between  these  two  solutions  to  derive  a  linear 
equation 

%N  f  a  uN  =  b 
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Fig.  5-10  Mesh  Configuration 


for  prODerties  behind  tne  shock;  (8)  solve  the  two  equations  from  step  (7)  ana 
step  (5)  for  the  shock  properties  at  SN.  The  accuracy  of  the  scheme  is 
furtner  enhanced  by  correcting  the  location  of  tne  shock  at  SN  based  on  an 
average  slope  between  S  and  SN  and  repeating  the  sequence  1-3  above. 

The  discontinuity  point  calculation  procedure  for  three  dimensional  flow 
requires  additional  analysis.  The  procedure  presented  for  two  dimensional 
flow  takes  full  advantage  of  supersonic  flow  theory.  "This  theory  provides  a 
direct  relationship  between  the  geometric  flow  angle#  and  the  Prandtl  Meyer 
angle  along  characteristic  lines.  In  three  dimensional  flow,  there  is  no 
simple  direct  relationship  between  flow  orientation  and  Prandtl -Meyer  angle. 
The  flow  direction  at  any  point  is  described  by  two  angles,  rand  o  for 
example,  and  is  not  unique.  In  the  past  investigators  (Refs.  5-6,  7,  8,9)  have 
employed  "reference  plane"  characteri sties  to  overcome  this  complication.  In 
this  method,  a  plane  is  prescribed  a  priori  which  is  oriented  in  a  direction 
convenient  to  the  calculation.  That  is,  the  direction  chosen  simplifies 
program  logic  and  reduces  the  need  for  interpolation  to  a  minimum.  The 
components  of  the  velocity  vector  in  these  planes  are  treated  as  two 
dimensional  flows  and  are  evaluated  based  on  a  method  of  characteri st i cs 
analysis  in  the  reference  plane.  Velocity  components  out  of  the  reference 
plane  are  computed  by  an  auxiliary  equation  -  usually  the  transverse  momentum 
equation.  Velocity  gradients  and  the  transverse  velocity  components  appear  in 
the  inhomogeneous  or  forcing  term  in  the  reference  plane  characteristic 
equations.  In  the  present  calculation  procedure,  a  new  method  is  described. 

It  is  based  on  a  theoretical  analysis  which  demonstrates  that  there  are  in 
fact  distinguished  planes,  in  a  local  sense,  in  which  the  flow  is  described  by 


equations  in  two  space  dimensions. 


In  three  dimensional  supersonic  flow,  the  zone  of  dependence  of  a  poi  n 


in  the  flow  is  the  upstream  Mach  cone.  Each  generator  of  this  cone  is  a 
oicharacteri Stic.  It  was  demonstrated  in  Ref.  5-10  that  two  of  the  infinite 
number  of  bicharacteristics  on  this  Mach  cone  play  exactly  the  same  role  as  do 
the  two  characteri sties  in  axisymmetric  flow.  These  bi characteri sties  are 
along  the  intersection  of  the  Mach  cone  and  the  osculating  plane  of  the 
streamline  Figure  5-11.  The  osculating  plane  contains  the  tangent  vector 
and  principal  normal  to  the  streamline.  The  normal  to  this  plane  is  in  the 
binormal  (or  b)  direction.  There  are  no  velocity  components  normal  to  the 
osculating  plane.  In  fact  the  equations  of  motion  normal  to  this  plane  reduce 
to  dp/db  =  0  (p  is  the  static  pressure  and  b  is  the  coordinate  in  the  bi normal 
direction).  In  the  osculating  plane  there  is  the  familiar  relationship 
between  changes  in  flow  angle  and  changes  in  Prandtl  Meyer  angle  that  exists 
for  axisymmetric  flow.  Only  the  gradient  in  the  b  direction  of  velocity 
component  in  the  b  direction  appears  in  the  equations.  The  entire 
computational  ^ rheme  devised  for  the  discontinuity  point  calculation  for  two 
dimensional  (or  axi symmetric )  flows  can  be  employed  in  the  three  dimensional 
flow  in  the  osculating  plane.  In  general  the  binormal  direction  is  given  by 
b  =  VP  x  q  (d  is  the  local  velocity  vector,  b  is  a  vector  parallel  to  the  unit 
vector  b).  This  formula  is  based  on  the  properties  that  the  b  vector 
be  normal  to  both  the  tangent  vector  and  acceleration  vector  of  a  fluid 
element,  q  is  tangent  to  the  streamline  andyp  is  parallel  to  the 
acceleration  vector.  Thus  the  orientation  of  the^  vector  or  alternatively 
the  orientation  of  the  osculatinq  plane  are  determined  alonq  with  the 
solution. 
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To  demonstrate  the  use  of  this  theory  the  three  dimensional  shock  point 
calculation  is  outlined  pe'ow.  Figure  5-12  shows  a  schematic  of  the  shock 
point  calculation.  A1  1  flow  properties  are  known  at  station  z  including  the 
shock  velocity  -n  the  /  direct ’on.  The  procedure  follows  exactly  the  steps 
outlined  in  the  tw«'  nmensional  calculation  described  in  the  previous 
oaragraohs  with  tne  following  additions.  The  shock  position  at  all  points  at 
z  --  Ac  are  computed  as  a  first  estimate  by  simple  forward  integration.  Thus 
YSN  =  VS  +  (dYS;dz'Az 

where  YSN  (YS)  is  tne  value  of  the  shock  height  at  2  +  Az  (z)  on  the  vertical 
grid  line  through  mesh  point  B  (N)  and  dYS/dz  is  the  shock  velocity  in  the  y 
direction  at  station  z  on  the  same  grid  line.  After  computing  the  new  location 
of  the  shock  on  each  mesh  line  the  shock  cross  flow  angle  can  be  determined  by 
a  simple  centered  difference  of  the  locations  at  SN  +  and  SN  ‘  (see  Fig. 

5-12).  The  only  unknown  parameter  for  the  shock  geometry  at  point  SN  is 
dYSN/dz  or  the  shock  speed  at  station  z  +  Az.  As  in  the  two  dimensional  case 
(step  6)  two  values  are  assumed  for  the  shock  speed.  In  the  3D  case  the  two 
velocity  vectors  associated  with  these  two  assumed  shock  normals  are  used  to 
define  the^  direction  (Fig.  5-13).  The  remainder  of  the  calculation 
procedure  is  identical  with  the  2D  procedure.  It  is  not  clear  that  tn’s  is 
only  method  for  computing  the  ^  vector.  The  method  is  linked  to  the  way  the 
shock  cross  angle  is  computed.  This  is  nonunique  and  the  influence  of 
variations  of  this  step  on  calculated  results  should  be  studied.  The  method 
presented  is  self  consistent  and  achieves  accurate  and  reliable  results  when 
combined  with  the  rest  of  the  computational  procedure. 
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Fig.  5-12.  Schematic  Diagram  Showing  Notation  for  Three-Dimensional 
Shock  Computation 
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A  computation  employing  the  three  dimensional  floating  discontinuity 
program  was  achieved  for  the  flow  field  created  by  the  impingement  of  two 
uniform  rectangular  jets.  This  geometry  was  chosen  so  as  to  minimize  the 
complexity  brought  about  by  nonuniformities  in  the  underexpanded  plume.  The 
geometry  describing  the  initial  calculation  is  shown  in  Fig.  5-14.  Two 
uniform  Mach  3.0  plumes  of  rectangular  cross  section  impinge  at  30°. 
Impingement  shocks  spread  across  the  plumes  (top  view)  to  make  the  two  flows 
parallel.  This  results  in  pressure  above  the  background,  and  the  flow  spreads 
laterally  (side  view)  to  relieve  this  overpressure.  The  cross-section  shown  in 
Fig.  5-14  is  characteristic  of  the  calculated  results.  The  impingement  shocks 
are  slightly  curved  and  are  bounded  by  the  free  jet  boundary.  The  pressure 
boundaries  spread  laterally  in  a  vee  shaped  pattern.  A  typical  cross  section 
from  the  calculation  is  shown  in  Fig.  5-15.  (Only  one  fourth  of  the  total 
cross  section  is  shown  because  the  flow  has  bilateral  symmetry).  The  flow 
from  the  undisturbed  plume  passes  downward  through  the  impingement  shock  and 
jumps  in  pressure.  The  impingement  shock  intersects  the  undisturbed  plume 
boundary  in  a  complex  interaction  involving  a  sonic  shock  condition  and  a 
sonic  shock  condition  and  a  centered  Prandtl  Meyer  fan  with  the  combined 
result  of  no  pressure  change  along  the  pressure  boundary.  Figure  5-16  is  a 
schematic  of  the  flow  near  that  point.  The  flow  is  decomposed  into  components 
tangential  and  normal  to  the  shock  wave/plume  boundary  intersection  line. 

(see  Section  6  of  this  paper  for  detailed  discussions  relating  to  this  type  of 
analysis).  In  the  plane  normal  to  the  intersection  line  the  shock  wave  is 
sonic  at  the  point  of  impingement.  Attached  to  the  impingement  point  is  a 
centered  expansion  which  reduces  the  pressure  back  to  ambient.  The 
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Fig.  5-14  Schematic  of  Geometry  for  the  Impingement 
of  Two  Rectangular  Plumes 
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Fig-  5-15  Typical  Cross  Section  at  2  =  2.12  Including  Isobars 
for  the  Sample  Calculation 


Fig.  5-16.  Details  of  Impingement  Shock/Plume  Boundary  Interac 
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fundamental  unknown  in  the  calculation  of  this  point  is  the  speed  of  the 
impingement  point.  The  calculation  procedure  is  similar  to  ordinary  shock 
points.  Flow  properties  behind  the  impingement  shock  are  matched  to  the 
interior  flow  by  a  characteristic  relationship  in  a  local  osculating  plane. 
This  procedure  is  an  approximation  and  must  be  examined  in  the  light  of  the 
more  accurate  flow  solutions  that  are  discussed  in  Section  6. 

The  isobars  for  the  cross  section  in  Fig.  5-15  show  that  the  flow  has 
nearly  the  undisturbed  two  dimensional  impingement  shock  value  at  the 
centerline  (see  Fig.  5-17).  The  decay  to  background  pressure  takes  place 
across  the  entire  flow  and  is  most  rapid  in  the  vicinity  of  the  shock/boundary 
intersection  point.  Figure  5-17  shows  the  symmetry  plane  pressure  profile  and 
the  cross  sectional  view  z  =  2.12  (the  plume  half  width  is  unity).  There  is  a 
region  of  near  constant  pressure  developing  at  the  outer  fringes  of  the 
pressure  boundary  as  would  be  expected.  There  are  some  "wiggles"  in  the 
pressure  in  this  zone  that  are  most  probably  due  to  the  low  number  of  mesh 
points  used  in  the  calculation  (Fig.  5-15  shows  the  exact  mesh  employed). 
Figure  5-18  (a)  shows  the  calculated  development  of  the  cross  sections  for  the 
impingement  region  as  a  function  of  distance  downstream  of  the  impingement 
line.  Each  profile  is  ten  calculation  steps  from  the  previous;  the  first 
being  at  Step  10.  The  pressure  boundary  develops  into  a  pointed  shape. 

Figure  5-18  (b)  shows  the  flow  profiles  superimposed  on  the  computational 
mesh.  The  impingement  shock,  pressure  boundary  and  interaction  cell  have 
crossed  many  mesh  points  and  have  attained  a  variety  of  conf l gurat ions  with  no 
apparent  breakdown  of  the  scheme.  An  interestino  comparison  is  made  in  Fig. 
5-19  with  calculations  reported  in  Ref.  5-11.  The  calculations  are  for  the 
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Fig.  5-17  Computed  Symmetry  Plane  Pressure  Profile  z  -  2.12 
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Fig.  5-19  Comparison  of  Present  Results  with  Scramjet  Plume  Calculation 


plume  boundary  of  a  scramjet  exhaust  employing  a  shock  capturing  technique. 

The  splitter  plate  produces  impingement  shocks  similar  in  geometry  to  those  in 
the  present  calculation.  The  comparison  which  is  meant  to  be  qualitative,  is 
quite  striking  in  that  even  the  irregular  shape  of  the  boundary  is  reproduced. 
This  gives  some  confidence  in  the  present  results,  however,  further  comparison 
with  other  three  dimensional  flow  calculations  would  be  useful. 
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SECTION  6 


THREE  DIMENSIONAL  PHENOMENA 

In  Section  3  the  flow  pattern  of  two  underexDanded  interacting  rocket 
plumes  was  discussed  and  the  general  features  of  the  possible  flow  fields  were 
pointed  out.  In  this  section  the  intersection  of  two  three  dimensional  shock 
surfaces  and  the  impingement  shock  liftoff  will  be  analyzed.  These  features 
have  no  counterparts  in  two  dimensional  supersonic  flow.  In  both  situations 
there  is  an  abrupt  change  in  configuration  because  attached  shock  solutions 
are  no  longer  viable.  The  conf i gurational  change  takes  place  through  a 
conical  solution  centered  at  the  transition  point.  Detailed  understandi ng  of 
the  nature  of  these  transitions  and  the  ability  to  compute  them  accurately  is 
required  in  the  subject  problem.  Many  other  three  dimensional  supersonic  flows 
of  current  interest  require  understanding  of  the  shock  intersection  problem 
and  other  three  dimensional  features  which  are  peculiar  to  those  problems. 
These  transition  problems  are  generally  complex.  The  shock  wave  pattern 
formea  by  two  intersecting  wedges  is  an  example  of  a  similar  conical  fiow  that 
is  less  complex  than  the  two  transitions  being  considered.  The  computation  of 
the  "corner  flow"  problem  has  been  the  subject  of  research  for  many  years  and 
is  the  subject  of  several  papers  (Refs  5-1  to  6-3)  and  was  recently  the 
subject  of  a  Ph.D.  thesis  (Ref.  6-4)). 

The  intersection  process  that  occurs  between  two  three  dimensional  shocx 
surfaces  will  be  discussed  first  .  A  single  three  dimensional  shock  surface 
has  a  certain  amount  of  arbi tra ri ness  in  its  description  that  does  not  exist 
in  the  two  dimensional  counterpart.  In  both  cases,  the  jump  conditions  for 


the  thermodynamic  state  variables  can  be  reduced  to  the  one  dimensional 
Rankine-Hugoniot  relationships  employing  the  Mach  number  normal  to  the  shock 
wave  (defined  by  velocity  component  normal  to  the  wave  and  the  undisturbed 
sound  speed).  The  tangential  component  of  the  oncoming  velocity  vector  is 
unchanged  by  the  shock  wave.  In  two  dimensions,  this  leads  to  an  unambiguous 
relationship  between  pressure  ratio  and  flow  deflection  across  a  shock  wave 
and  the  definition  of  a  shock  polar.  In  three  dimensions,  it  is  possible  to 
arbitrarily  decompose  the  tangential  velocity  vector  into  two  components  in 
the  surface  of  the  wave.  One  of  these  components  can  be  combined  with  the 
normal  velocity  component  to  form  a  velocity  vector  oblique  to  the  wave.  The 
jump  condition,  for  the  thermodynamic  variables,  can  be  determined  by  the 
corresponding  two  dimensional  oblique  shock  relationships.  The  velocity 
downstream  of  the  shock  surface  is  found  by  adding  to  the  oblique  shock 
solution  the  remaining  component  of  tangential  velocity  (which  is  unchanged  by 
the  wave).  Because  the  decomposition  of  the  tangential  velocity  vector  is 
arbitrary,  a  three  dimensional  shock  wave  surface  can  be  described  by  an 
infinity  of  two  dimensional  oblique  shocks  with  an  additional  tangential 
velocity  component.  For  a  given  surface  normal  vector  to  the  shock  the 
downstream  result  is  always  identical. 

The  intersection  of  two  shock  waves  in  three  dimensional  flow  takes  place 
along  an  intersection  line  whereas  the  two-dimensional  counterpart  takes  place 
at  a  point.  The  properties  of  the  interaction  depend  not  only  on  the  shock 
strengths  but  also  the  local  orientation  of  the  intersection  line.  When  the 
flow  is  supersonic  downstream  of  the  interaction  the  intersection  line  is  at 
the  leading  edge  of  the  transmitted  and  reflected  waves.  The  transmitted  and 
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reflected  waves  originate  at  the  intersection  line  in  the  same  manner  the 
attached  shock  surfaces  originate  at  the  leadinq  edges  of  swept  wings  in 
supersonic  flow.  Thus  the  sweep  of  the  intersection  line  plays  a  strong  role 
in  the  interaction.  Because  a  three  dimensional  shock  surface  can  cause  flow 
deflection  in  any  azimuthal  direction  it  cannot  be  simply  described  as  a  first 
(upward  deflection)  or  a  second  (downward  deflection)  “family"  shock  as  a  two 
dimensional  shock  can.  In  two  dimensional  flows  the  intersection  of  two  shock 
waves  can  be  categorized  according  to  whether  both  shocks  are  of  the  same  or 
opposite  family.  In  either  case  there  are  a  variety  of  shock  configurations 
that  can  result  and  the  subject  is  an  established  segment  of  supersonic  flow 
theory  (see,  for  example,  Refs.  6-5,6).  While  a  three  dimensional  shock 
cannot  be  described  simply  as  belonging  to  one  of  two  families,  the 
interaction  process  of  two  shock  surfaces  can  still  be  categorized  in  this 
manner.  At  a  point  along  the  intersection  line  each  of  the  two  waves  can 
rotate  the  flow  in  the  same  direction  or  in  opposite  directions  relative  to  an 
imaginary  axle  parallel  to  the  intersection  line.  Thus,  the  simple 
classification  employed  in  the  two  dimensional  theory  can  be  employed  locally 
for  three  dimensional  shock  surfaces. 

In  the  following  discussion  the  evolution  of  a  three  dimensional 
shock/shock  interaction  process  is  elucidated  by  considering  the  model  problem 
of  a  circular  conical  shock  intersecting  a  planar  shock  surface.  This  model 
problem  has  the  highly  desirable  feature  of  constant  shock  strengths  (in  terms 
of  pressure  ratios)  at  all  points  along  the  geometric  intersection  line  which 
is  a  true  hyperbola.  The  case  considered  here  has  a  regular  reflection 
process  at  the  leading  edge  of  the  intersection  line.  When  viewed  along  the 
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intersection  line,  the  initial  interaction  is  locally  that  of  two  same  family 
shock  waves  intersecting  and  there  is  a  simple  transmitted  and  reflected  wave 
pattern.  Downstream  of  the  leading  edge  the  sweep  angle  of  the  intersection 
line  increases,  giving  rise  to  decreasing  apparent  normal  Mach  numbers  and 
increasing  apparent  deflection  angles  so  that  at  some  point  the  simple  regular 
reflection  process  is  no  longer  possible.  A  detailed  discussion  of  the 
reasons  for  this  is  given  employing  the  shock  polar  diagrams  in  the 
pressure/deflection  plane.  A  Mach  (or  irregular)  reflection  process  develops 
in  which  an  additional  shock  segment  bridges  the  span  between  the  two  incoming 
waves.  This  shock  segment  moves  ahead  of  the  geometrical  intersection  of  the 
two  incoming  waves  and  there  are  now  two  shock  triple  points  characterizi ng 
the  intersection.  The  cross  flow  behind  this  additional  shock  segment  is 
subsonic  relative  to  the  triple  points  and  is  analogous  to  the  Mach  reflection 
that  develops  in  the  three  dimensional  compression  corner  flow  field.  The 
propagation  of  an  irregular  reflection  pattern  thus  involves  the  intersection 
of  three  shock  surfaces  -  the  two  incident  waves  and  the  Mach  reflection  wave 
which  is  generated  by  the  intersection. 

The  complex  nature  of  the  intersection  of  shock  wave  surfaces  in  three 
dimensional  flows  can  be  illuminated  by  considering  the  model  flow 
configuration  shown  in  Fig.  6-1.  A  thin  wedge  in  a  supersonic  stream  produces 
a  supersonic  oblique  planar  shock  surface  and  a  uniform  flow  parallel  to  the 
wedge.  A  cone  is  aligned  with  this  flow  giving  rise  to  a  conical  shock 
surface.  The  conical  shock  and  the  planar  shock  intersect  along  a  hyperbolic 
line  in  the  plane  of  the  oblique  shock  (Fig.  6-1 ( b ) ) .  In  cross-sectional 
view,  the  circular  conical  shock  intersects  the  planar  oblique  shock  at  two 
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bl  INTERSECTION  LINE  IN  PLANE  OF  THE 
WEDGE  SHOCK 


c)  CROSS  SECTIONAL  VIEW  OF  SHOCK  SYSTEM 
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Fig.  6-1  Intersection  of  a  Planar  Wedge  Shock  and  a  Circular 
Conical  Shock 
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points.  The  wave  system  that  joins  these  two  points  can  be  characterized  as 
either  a  simple  transmitted  and  reflected  wave  system  (regular)  or  a  more 
complex  system  containing  the  equivalent  of  a  Mach  reflection  (irregular). 
These  are  shown  on  the  right  hand  side  and  left  hand  side  of  Fig.  6-1 ( c ) 
respectively.  The  pattern  that  prevails  (regular  or  irregular)  depends  on  the 
location  along  the  intersection.  Initially,  the  sweep  angle  is  zero  and  the 
interaction  is  regular  (for  the  cases  considered  here  with  supersonic 
downstream  velocity).  Downstream,  as  the  sweep  angle  increases,  a  transition 
occurs  and  the  shock  pattern  becomes  irregular.  In  the  regular  reflection  a 
local  two  dimensional  flow  analysis  centered  on  the  intersection  line  can  be 
employed  to  determine  the  local  solution.  There  is  a  simple  transmitted  and 
reflected  wave  and  contact  surface  that  leaves  the  point  of  intersection 
(right-hand  side  of  Fig.  6-1 ( c ) -  The  cross  flow  relative  to  the  intersection 
line  behind  the  reflected/transmitted  waves  is  supersonic.  Therefore,  the 
intersection  line  can  be  viewed  as  the  leading  edge  of  these  shocks  and  the 
solution  is  independent  of  downstream  properties  .  In  the  Mach  (irregular) 
reflection  case,  the  resulting  shock  system  is  actually  ahead  of  the 
geometrical  intersection  of  the  cone  and  wedge  shocks  (left  hand  side  of  Fig. 
6-1  ( c ) .  There  are  now  two  shock  triple  points,  and  the  cross  flow  is,  in 
general,  subsonic;  thus,  flow  properties  from  the  high  pressure  side  determine 
the  propagation  velocity. 

In  the  case  of  the  regular  reflection  (shown  schematically  in  Fig.  6-2, 
the  solution  is  arrived  at  by  decomposing  the  free  stream  velocity  into  the 
components  normal  and  tangential  to  the  intersection  line.  The  tangential 
component  of  velocity  is  parallel  to  all  the  wave  surfaces  by  virtue  of  the 
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Fig.  6-2  Isometric  View  of  Regular  Shock  Pattern  for 
the  Wedge-Cone  Shock  System 
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geometry  of  the  intersection.  Figure  6-3  details  a  point  along  the 
intersection  line  showing  the  two  incoming  waves  and  the  transmitted  and 
reflected  wave  surfaces.  In  the  plane  containing  the  normal  component  of 
velocity  (Fig.  6-4),  the  flow  pattern  is  identical  to  that  of  the 
two-dimensional  flow  containing  the  intersection  of  two  shock  waves  of  the 
same  family.  The  incident  shock  waves  are  known  in  strength  and  do  not 
interact  until  they  intersect  (the  flow  behind  them  is  supersonic).  The 
solution  for  the  transmitted  and  reflected  wave  strength  is  achieved  by 
determining  the  deflection  of  the  contact  surface  63  (Fig.  6-4)  to  match  the 
pressures  in  regions  (3)  and  (4)  behind  the  transmitted  and  reflected  waves, 
respectively.  This  can  be  achieved  by  a  simple  iteration  procedure.  The 
process  is  best  understood  by  studying  the  flow  process  in  the  pressure/flow 
deflection  hodograph  plane  (Fig.  6-5).  Each  point  on  curve  I,  the  shock  polar 
associated  with  the  upstream  velocity  (the  velocity  component  normal  to  the 
intersection  line  in  the  three-dimensional  case)  corresponds  to  conditions  on 
the  high  pressure  side  of  the  shock  wave.  The  entire  polar  corresponds  to  all 
shock  wave  angles  ranging  from  the  local  Mach  angle  to  a  normal  shock  (for 
both  positive  and  negative  deflections).  The  pressure  in  zone  (1)  is 
determined  simply  by  moving  to  the  defl ection  d j ,  on  the  abscissa.  For  the 
Mach  number  of  the  flow,  in  zone  (1)  a  shock  polar,  curve  II,  can  be 
constructed  at  point  (1),  and  the  pressure  at  (2)  is  determined  corresponding 
to  the  deflection  $2-  The  shock  polar  at  point  (2),  curve  III  completes 
the  solution.  Polars  I  and  III  intersect  at  a  point  which  is  denoted  as  point 
(3)  along  polar  I  and  (4)  along  III.  Points  (3)  and  (4)  are  coincident  in  the 
pressure/deflection  plane  but  represent  different  velocities  and  entropy  (or 
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Fig.  6-3  Details  of  Regular  Intersection  Pattern  of  the  Wedge 
Shock  and  Cone  Shock  Along  the  Intersection  Line 


Fig.  6-4  Regular  Shock  Reflection  Pattern  in  the 
Plane  Normal  to  Intersection  Line 
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Fig.  6-5  Pressure  Deflection  Plane  Showing  Shock 
Polar  Configuration  for  Regular  Reflection 
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density).  3y  the  nature  of  the  construction,  point  (3/4)  is  at  the  desired 
condition  which  brings  the  flow  behind  the  transmitted  and  reflected  waves  to 
identical  pressures  and  flow  deflections.  The  regular  wave  system  depicted  in 
Fig.  6-5  has  a  weak  compression  wave  (the  arc  (2)  -  (3))  as  the  reflected 
wave. 

The  reflected  wave  between  polars  I  and  II  can  be  an  expansion  or 
compression  wave  depending  on  the  orientation  of  polars  I  and  II.  There  are 
four  possible  intersection  patterns  of  polars  I  and  II  (Refs.  6-5,  6-6)  shown 
in  Fig.  6-6  that  are  characterized  by  the  number  of  points  the  two  curves 
intersect.  Whenever  point  (2)  on  polar  II  (see  Fig.  6-5)  is  at  a  lower 
pressure  than  a  point  on  curve  I  for  the  same  deflection,  the  reflected  wave 
is  a  compression  (shock)  as  described  in  the  preceeding  paragraphs.  However, 
in  the  m  -  0,  m  =  2  and  a  portion  of  the  m  -  3  case  (m  is  the  number  of 
intersection  points  of  the  two  curves),  the  geometry  is  reversed  and  the 
solution  for  the  reflected  wave  is  a  simple  centered  expansion,  as  shown  in 
Fig.  6-7.  It  can  be  shown  (Refs.  6-5,  6-6)  that  the  intersection  pattern  is 
dependent  on  free  stream  Mach  number  and  deflection  angle  6^,  (as  well  as 
the  ratio  of  specific  heats  for  the  gas).  In  any  of  these  cases  the  wave 
pattern  becomes  irregular  when  the  shock  polar  II)  does  not  intersect  polar  I 
(or  in  the  case  of  a  reflected  exnansion  wave  when  the  characteri stic  curve 
does  not  intersect  I).  The  simple  all  supersonic  (cross)  flow  pattern 
associated  with  the  regular  reflection  is  no  longer  possible  and  a  more 
complex  pattern  emerges. 

In  the  model  flow  problem  for  specified  cone  and  wedge  angles  and  free 
stream  Mach  number  (the  conic a'l /wedge  shock  intersection  problem),  the  local 
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.  6-6  The  Four  Possible  Configurations  tor  Two 
Shock  Intersections 
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Pig,  6-7  Regular  Shock  Intersection  with  a  Reflected 
Expansion  Wave 


solution  depends  solely  on  the  sweep  angle  (see  Fig.  6-1)  of  the  intersection 

line.  This  is  because  the  pressure  ratios  across  each  shock  wave  are  constant 

and  independent  of  location.  In  addition  the  apparent  deflections  through  the 

two  incoming  waves  as  well  as  the  governing  Mach  number  normal  to  the 

intersection  line  are  all  functions  of  the  sweep  angle  because  the  flow  is 

uniform  upstream  of  the  wedge  shock.  Evolution  of  the  shock  system  transition 

from  regular  to  irregular  is  illustrated  in  Fig.  6-8.  These  figures  detail 

the  shock  polars  for  the  following  conditions:  =  3.0,  y  -  1.4,  wedge  angle 

10°,  wedge  pressure  ratio  Pw/Poo  =  2.054,  cone  angle  9C  3  20°,  cone 

pressure  ratio  Pc/Pw  =  1.95  or  Pc/Poo  =  4.005,  (see  Fig.  6-1  for 

definitions).  The  asymptotic  sweep  angle  of  the  intersection  line  (^  )  is 

60°.  Figure  6-8a  is  the  pressure/deflection  plane  for  the  initial 

intersection  of  the  shocks  (0=0).  In  this  case  the  polars  are  very  close 

! 

together  and  the  reflected  wave  is  a  very  weak  expansion  wave  which  cannot  be 
seen  on  the  scale  of  the  figure.  Downstream  along  the  intersection  line  Fig. 
6-8b  shows  the  shock  polar  pattern  for  0  =  52°.  Note  that  while  the  polar 
patterns  are  considerably  different,  the  pressure  at  points  (2)  and  (3)  which 
are  the  wedge  and  cone  shock  pressures  are  the  same  as  the  0=0°  case.  At 
this  sweepback  angle,  0=  52°,  the  normal  Mach  number  to  the  intersection  line 
has  decreased  to  2.15  from  the  free  stream  value  of  3.0.  The  apparent  wedge 
and  cone  shock  deflection  angles  which  were  10°  and  20.8°  initially  are  now 
13°  and  26.2°,  a  geometric  result  of  projecting  the  flow  fields  onto  the  plane 
normal  to  the  intersection  line.  The  decreased  normal  Mach  number  and 
increased  apparent  deflection  moves  the  solution  toward  the  upper  right  hand 
corner  of  polar  I.  A  closeup  of  this  upper  right  hand  portion  of  Fig.  6  -8b 
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Fig-  6-8  Pressure  Deflection  Diagram  for  the  Intersection  of  a 
Conical  Shock  and  a  Wedge  Shock  (Sheet  1  of  2) 


is  shown  in  Fig.  6-8c.  This  case  is  interesting  because  point  (4),  the  flow 
behind  the  transmitted  wave,  is  on  the  subsonic  portion  of  polar  I.  The  cross 
flow  becoming  subsonic  is  a  signal  that  there  is  an  incipient  change  in 
structure.  In  fact,  a  small  distance  downstream  where  0  =  55°  (Fig.  6-8d) 
polar  III  no  longer  intersects  polar  I,  and  the  regular  four  wave  pattern  that 
had  prevailed  is  now  not  possible.  A  closeup  of  the  interesting  portion  of 
this  figure  is  shown  in  Fig.  6-8e* 

The  subsequent  flow  pattern  is  characterized  by  a  five-wave  intersection 
made  up  of  two  three  wave  intersections,  as  shown  in  Fig.  6-9a.  This  shock 
pattern  is  qualitatively  similar  to  the  one  that  prevails  in  the  internal 
compression  corner  flow  field  (Refs  6-1 ,2,3,4).  The  major  qualitative 
difference  is  that  the  two  contact  surfaces  forming  the  triangular  zone  behind 
the  Mach  reflection  shock  meet  at  a  point,  denoted  b  on  Figure  6 -9a,  on  a 
contact  surface.  In  the  compression  corner  flow  field  this  point  is  in  the 
corner.  It  is  not  possible  to  characterize  the  entire  interaction  by  a  single 
shock  polar  pattern  because  there  are  now  two  intersection  lines.  However,  in 
order  to  visualize  the  process  the  velocity  field  normal  to  the  interaction  is 
approximated  as  constant  and  uniform.  In  Fig.  6-9b  the  flow  in  the  pressure 
/deflection  plane  would  then  be  as  follows  (Ref.  6-6)  describes  the  two 
dimensional  counterpart):  (a)  points  (1)  and  (2)  are  at  the  wedge  and  cone 
shock  pressures  as  before;  (b)  shock  polars  I  and  II  intersect  in  only  one 
place  (for  this  case)  and  that  determines  the  pressure  and  deflection  at  the 


*  The  change  in  structure  might  come  at  the  point  where  the  cross  flow  becomes 
sonic.  It  is  not  clear  at  present  exactly  where  the  transition  occurs. 
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SUM* 


Intersection  line  labeled  A;  (c)  a  similar  process  takes  place  at  intersection 
B  -  the  junction  of  zones  (1),  (2),  (3)  and  (5)  -  where  B  (in  Fig.  6-9b) 
represents  the  only  possible  shock  transition  starting  from  conditions  (1)  and 
(2);  and  (d)  where  the  segment  of  the  shock  wave  between  A  and  B  is  a  curved 
shock  (shown  as  a  bold  line  in  Fig.  6-9b). 

A  schematic  of  the  expected  transition  pattern  is  shown  in  Fig.  6-10. 
Relative  to  the  transition  or  breakup  point  the  flow  is  conical  in  nature. 

The  cross  flow  is  subsonic  in  the  zone  immediately  behind  the  shock 
interaction  and  is  thus  governed  by  elliptic  partial  differential  equations. 
Thus  the  propagation  speeds  of  the  two  triple  points  leaving  the  transition 
zone  can  only  be  determined  as  part  of  the  solution  of  the  entire  transition 
zone.  An  interesting  and  important  aspect  of  the  flow  pattern  is  that  the 
pressure  levels  prevailing  behind  the  interaction  shocks  are  higher  than  the 
levels  attained  by  the  regular  reflect  pattern  just  before  it  breaks  apart 
(shown  as  the  X  on  polar  I  in  Fig.  6-9b).  Thus,  at  the  point  along  the 
intersection  line  where  the  shock  pattern  transitions  from  regular  (four-wave) 
to  irregular  (five-wave)  the  flow  pattern  must  include  a  three  dimensional 
expansion  zone  immediately  behind  the  shock  wave  to  equilibrate  the  pressure. 
Again  the  strength  and  distribution  of  this  expansion  zone  are  determined  as 
part  of  the  complete  solution  of  the  transition  zone. 

The (impingement)  shock  lift  off  problem  is  another  three  dimensional 
transition  problem  that  occurs  in  the  multi  nozzle  flow  field.  Between 
stations  denoted  by  0  and  1  Figure  3-2(a)  is  the  initial  portion  of  the 
impingement  flow  field  and  the  origin  of  the  impingement  shock  surface.  The 
flow  pattern  in  this  region  can  be  studied  by  examining  the  simpler  flow  in 
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which  two  uniform  circular  jets  impinge.  Figure  6-11  shows  the  overall  flow 
configuration  (leaving  out  the  interaction  for  clarity).  Two  circular  jets 
with  uniform  properties  impinge  at  the  symmetry  plane.  The  intersection  of 
the  plume  boundary  and  the  symmetry  plane  is  an  ellipse.  The  shock  pattern 
develops  as  is  shown  in  an  isometric  view  in  Fig.  6-12.  In  the  cases  that  are 
considered  here  the  shock  wave  is  initially  attached  at  the  origin.  The 
initial  shock  pattern  when  viewed  as  cross  sections  have  the  appearance  of  arc 
segments  anchored  at  the  ends  to  the  intersection  line.  At  some  point  along 
the  intersection  line  this  pattern  is  no  longer  possible  and  the  shock  ends 
lift  off. 

The  shock  pattern  is  shown  in  Figure  6-13  along  with  the  details  of  the 
flow  at  the  shock  leading  edge.  In  a  manner  completely  analogous  to  that 
employed  to  determine  the  shock  wave  strength  at  the  leading  edge  of  a  swept 
wing  in  supersonic  flow  (c.f.  Ref.  6-7)  the  shock  is  examined  in  a  plane 
perpendicular  to  the  intersection  line.  Again  the  component  of  velocity 
tangential  to  the  intersection  line  is  unaltered  by  the  shock  surface.  These 
sections  appear  adjacent  to  each  cross  section  in  Fig.  6-13  .  Near  the 
leading  edge  the  Mach  number  normal  to  the  intersection  is  large  enough  and 
the  deflection  ( 6  ^)  small  enough  so  that  there  is  an  attached  shock 
solution  as  shown  by  the  adjacent  pressure  deflection  shock  polar  diagram. 
Downstream  at  station  2  the  normal  Mach  number  (Mn2)  has  decreased  and  the 
apparent  deflection  5 2  has  increased.  The  shock  polar  has  diminished  in 
size  and  the  required  deflection  (62)  is  getting  near  the  maximum  deflection 
possible  for  Mn2*  At  station  3  the  shock  has  already  lifted  off  the 
symmetry  plane.  This  is  required  because  at  Mn3  the  deflection  cannot 
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Fig.  6-11  Schematic  Diagram  Showing  the  Intersection 
of  Two  Circular  Jets 


6-23 


3  0 


CROSS  SECTION  OF  SHOCK  WAVE  AT  LEADING  PRESSURE/DEFLECf  ON 

SHOCK  PATTERN  EDGE  VI EWED  ALONG  IMAGE  OF  FLOW 

INTERSECTION  LINE 

Fig.  6-13  Cross  Section  Details  of  Impingement  Shock  and  Corresponding 
Hodograph  Diagrams 


be  achieved  by  an  attached  shock.  Note  that  Mn3  and  63  are  determined 
relative  to  the  geometric  intersection  line  of  the  plume  and  the  symmetry 
plane.  The  actual  leading  edge  of  the  shock  wave  moves  up  along  the 
undisturbed  plume  boundary.  The  flow  perpendicular  to  this  leading  edge  or 
intersection  line  must  be  at  least  sonic  (see  Section  3)  because  the  velocity 
must  be  supersonic  to  reexpand  and  match  the  pressure  on  the  plume  boundary. 

In  the  neighborhood  of  the  liftoff  point  the  flow  pattern  is  conical. 

With  the  liftoff  point  as  the  origin  of  a  spherical  coordinate  system  the  flow 
pattern  develops  along  spherical  rays  and  is  geometrically  similar  when  scaled 
to  the  distance  from  the  liftoff  point.  A  schematic  of  this  conical  flow  is 
shown  in  Fig.  6-14.  The  orientation  and  relative  position  of  this  pattern  is 
shown  in  Fig.  6-12.  The  Mach  cone  that  leaves  the  liftoff  point  limits  the 
extend  of  its  influence.  The  flow  outside  this  cone  is  uneffected  by  the 
liftoff.  Inside  the  region  bounded  by  the  Mach  cone  the  cross  flow  (conical) 
is  subsonic  and  the  governing  conical  equations  are  elliptic.  Therefore  the 
solution  of  this  domain  rests  on  an  iterative  or  relaxation  procedure.  The 
segment  of  the  shock  between  AB  (Figs.  6-14  and  6-12)  has  the  following 
properties.  At  point  B  the  shock  rotates  the  flow  opposite  in  sense  from  that 
at  point  A.  At  point  A  and  at  points  1  and  2  upstream  (see  Fig.  6-13)  the 
shock  rotates  the  oncoming  flow  through  a  clockwise  deflection.  At  point  B 
the  deflection  is  counterclockwise.  Between  A  and  B  the  shock  goes  through 
the  spectrum  and  produces  no  rotation  (is  normal  in  the  conical  sense)  at 
point  C.  The  shock  produces  sonic  velocity  at  point  B  where  there  is 
anchored  a  centered  expansion  to  reestablish  the  background  pressure.  This 
expansion  produces  supersonic  cross  flow  so  that  a  portion  of  the  domain  shown 
in  Fig.  6-14  does  not  have  subsonic  cross  flow. 
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SECTION  7 


A  NEW  METHOD  FOR  TRACKING  DISCONTINUITIES  IN  NUMERICAL  CALCULATIONS 


The  multinozzle  plume  flow  field  is  an  example  of  a  complex  three 
dimensional  supersonic  flow  field.  These  flow  fields  by  their  very  nature 
have  solutions  which  are  discontinuous  or  have  discontinuous  derivatives. 

This  is  reflected  in  the  mathematical  description  of  the  flow  field  by  a  set 
of  hyperbolic  partial  differential  equations.  Equations  of  this  type  have 
wave  like  solutions  which  permit  the  propagation  of  functions  which  have 
discontinuous  derivatives.  In  addition  these  equations  admit  solutions  which 
are  discontinuous.  For  the  Euler  equations  the  jump  conditions  along  these 
discontinuities  are  the  familiar  Rankine-Hugoniot  conditions.  Two  regions  of 
continuous  solution  can  be  joined  by  a  jump  along  a  shock.  A  numerical 
solution  procedure  for  supersonic  flow  problems  should,  therefore,  include 
within  it  the  ability  to  deal  with  these  properties.  Section  4  of  this  paper 
describes  a  "floating  fitted  discontinuity"  scheme  that  tracks  each 
discontinuity  in  detail,  and  computes  the  jump  properties  exactly.  This 
method  as  outlined  requires  a  large  amount  of  program  logic  to  handle  complex 
geometries.  In  this  section  a  new  approach  to  tracking  di sconti nuites  is 
examined.  The  algorithm  employs  the  normal  functional  approximations  (Taylor 
series)  at  all  regular  mesh  points.  At  mesh  points  which  are  recognized  as 
being  adjacent  to  discontinuities  the  appropriate  flow  properties  are 
approximated  by  discontinuous  functions. 

There  is  an  underlying  assumption  in  almost  all  finite  difference 
algorithms  -  the  unknown  function  can  be  approximated  locally  by  a  Taylor 
series  expansion.  In  fact,  finite  difference  schemes  are  categorized 
primarily  on  this  basis.  A  scheme  is  said  to  be  second  order,  for  example, 
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if  the  difference  solution  and  the  exact  solution  expanded  in  a  Taylor  series 
match  to  second  order.  Supersonic  flows  solution  generally  cannot  be 
approximated  by  Taylor  series  in  all  regions.  It  seems  reasonable  then  to  try 
to  employ  other  more  general  functions  in  portions  of  the  flow  where  there  are 
discontinuities.  As  an  example  the  Euler  equations  for  2D  flow  can  be  written 
in  vector  form*  as 

Ux  +  Vy  =  0  (7-1) 

At  a  shock  wave  the  values  of  U  upstream  of  the  shock  (U')  and  downstream  of 
the  shock  (U+)  are  related  by 

[U]  -  [V]/W  =  0  (7-2) 

where  [  ]  means  jump  across  the  shock  and  W  is  the  shock  propagation  speed. 

The  exact  solution  in  the  vicinity  of  a  shock  wave  through  the  origin  is 
U  -  U"  +  [U]  H  (x  -  y/W) 

V  ■  V“  -  W[U]  H  (x  -  y/W)  (7-3) 

where  H(z)  is  the  unit  step  function  given  by 
H  (z)  »  0  z  <  0 

=1  z  >  0  (7-4) 

Simply  substituting  (7-3)  shows  that  it  is  a  solution  (7-1)  (in  the  context  of 
the  existence  of  the  delta  function  fl  (z)  )  This  analysis  strongly  suggests 
that  solutions  for  U  in  the  vicinity  of  the  shock  should  take  the  form 
U  =  IT  +  [U]  h  (z)  (7-5) 

where  h  (z)  is  a  chosen  function  that  has  the  desired  properties  of  H(z)  and 
can  easily  be  employed  in  a  numerical  scheme. 

-  /Pu  2\  /pv  \ 

*  for  the  Euler  equations  U=  [  p+  Pu  |  y=  (  2) 

I  puv  \  lP+pvc/ 

\puhQ  /  \pvhQ  / 
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Before  proceeding  it  is  interesting  to  note  some  of  the  properties  of  the 
"exact"  solution  for  a  shock  wave  Equation  (7-3).  Neither  U  or  V  are 
continuous  across  a  shock  wave.  U  +  WV  is  continuous 

U”  +  WV"  *  U+  +  WV+  (7-6) 

The  vector  V  is  a  function  of  U  only.  It  can  be  expressed  symbolically  as 
V  =  F(U)  so  that  V*  =  F(U+)  and  V"  =  F ( U~ ) .  Equation  7-6  is  then  a 
relationship  between  the  vector  U‘  and  U+ 

U-  +■  WF(IT)  *  U+  +  WF(U+)  (7-7) 

No  value  of  U  between  U"  and  1J+  can  be  substituted  on  the  right  hand  side 
of  Equation  (7-7).  In  shock  capturing  schemes  where  the  solution  vector  1) 
varies  smoothly  (or  not  so  smoothly)  between  U"  and  U+  in  the  reqion  of 
the  shock  jump  Equation  (7-7)  is  not  satisfied.  The  portion  of  the  flow  field 
between  U"  and  l)*  is  an  artifact  of  the  calculation  and  has  the  inherent 
error  associated  therein. 

As  a  test  bed  for  these  ideas  the  solution  of  a  supersonic  flow  field 
containing  a  contact  discontinuity  was  employed.  This  flow  field  was  chosen 
because  :  (1)  contact  discontinuities  spread  out  over  more  mesh  points  than 
shocks  and  (2)  it  represents  a  jump  in  entropy  (total  pressure)  only  - 
pressure  and  deflection  are  continuous  across  it.  In  the  following  paragraph 
the  use  of  a  ramp  function  to  model  a  jump  in  entropy  is  described  in  detail. 
Then  the  brief  discussion  of  the  numerical  calculation  procedure  is  given.  At 
mesh  points  away  from  the  jump  a  standard  first  order  finite  difference  scheme 
is  employed.  Results  of  two  test  computations  are  discussed  in  detail.  (In 
the  computer  code  stagnation  pressure  is  used  instead  of  entropy.) 


The  ramp  function  shown  in  Fig.  7-1  was  employed  in  order  to  implement 


some  of  these  concepts. 

II  z  <  0 

1  -  z/h  0  <  z  <  h 

0  h  <  z 

Figure  7-2  shows  a  portion  of  a  finite  difference  mesh  near  the  vicinity  of  a 

discontinuity  in  entropy  (or  stagnation  pressure).  In  the  mesh  inverval 
(j-1,  j+1 )  there  is  a  discontinuity  in  the  function  S.  If  the  functional  form 
of  S  in  this  mesh  interval  is  given  by  the  ramp  function  then 
S  -  Sj+1  +  (Sj_i  -  Sj+i)  R  (x  -  Xq j )  (7-9) 

where  x0j  =  xj  -  h (Sj_i~Sj )/(Sj_i  -Sj+i)  (7-10) 

or  A  =  (Sj.i  -  Sj)/(Sj_1  -  Sj+1) 
xoj  =  xj  h 

A  possible  interpretation  of  the  distribution  of  S  in  the  interval  is  shown 
as  a  dashed  line  in  Figure  7-2.  S  can  be  envisioned  to  have  a  step 
discontinuity  at  the  mid  point  of  the  ramp  function.  The  value  of  Sj  serves 
to  locate  this  point  between  j-1  and  j+1.  The  governing  equation  for  S  in  one 
dimensional  unsteady  flow  is* 

St  +  uSx  =  0  (7-11) 

Using  the  ramp  solution  of  equation  7-11  given  by  equation  (7-9)  the  value  of 
S  at  time  At  latter  than  is  shown  in  figure  1-2  is  given  by 


*  Discontinuous  solutions  of  this  equation  follow  the  same  form  as  equation 
7-3.  S  =  S'  +  [S]H(x-ut ) 
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SNj  =  Sj+i  +  (Sj_i  -  Sj+i)  R  (xj  -  u j  At  -  x0j)  (7-12) 

if  u  >  0 

SN j  +  ]_  =  +  +  (S j_p  Sj  +  ]_)  R  ( x j  -  iijAt  +  h  -  Xgj  )  (7-13) 

if  u  <  0 

SNj _ i  *  Sj  +  i  +  (Sj_!  -  Sj  +  i)  R  ( x j  -Uj At  -  h  -  Xg j )  (7-14) 

A  research  code  (KAYT1)  was  written  to  investigate  the  utility  of  the 
ramp  function  in  certain  problems.  The-  program  is  for  two  dimensional 
supersonic  flows  and  employs  a  finite  difference  calculation  procedure  using 
windward  differences  and  characteristic  form  for  the  equations.  The  purpose 
of  this  section  is  to  describe  the  results  using  the  ramp  function  so  only  a 
brief  outline  of  the  numerical  scheme  will  be  given.  The  primary  unknowns  are 
pressure,  flow  deflection  and  total  pressure.  The  governing  equations  are 


*x  +  Xj  0y  +  «W)(PX  +  Py)  =  0 

+a2  ey  -  (0/?m2)(px  +  a2  Py)  -  0 

Pox  +  *3  Poy-  0 

where  P  =  4n(p/pr),A1  =  tan(0+4),A2  -  tan  (0-4) , 
A3  =  tanfl,  and  PQ  =  ?n(PQ/Pr)  and  Pr  is 


(7-15) 

(7-16) 
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some  reference  value  of  pressure.  The  marching  direction  is  the  x  direction. 
Transverse  (y)  derivatives  in  each  equation  are  evaluated  by  either  forward  or 
backward  two  point  formulas  depending  on  the  sign  of  A,  multiplying  that 
term.  For  example,  in  equation  (7-15)  the  term  A]  Py  is  evaluated  at  mesh 
point  J  (where  y  =  h(J-l)  and  h  is  the  mesh  spacing)  as  follows 

(A./h)  (P(J)  -  P(J-l))  A.  >0 

A , P  »  1 

1  y  (A  /h)  (P(J+1)  -  P(J))  Xj<0 
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In  this  formula  all  values  of  the  pressure  are  at  the  known  station  x.  The 
derivatives  in  the  marching  direction  (x)  are  evaluated  by  two  point  forward 
derivatives  as  follows 

Px  =  (PN(J)  -  (P(J)  )/A x  (7-19) 

where  PN(J)  denotes  the  value  of  P  at  mesh  point  J  and  x  =  Ax.  Formulas 
(7-19)  and  (7-18)  are  used  for  all  derivatives  in  equations  (7-15),  (7-16)  and 

(7-17)  to  derive,  at  every  mesh  point,  two  equations  in  the  two  unknowns  PN 

and£N  and  an  equation  for  P0N.  The  code  overrides  the  basic  equation  for 
P0N  if  it  is  determined  that  there  is  a  discontinuity  in  P0  at  some  mesh 
point.  Then  the  determination  of  P0N  employs  the  ramp  function  algorithm  as 
outlined  in  the  previous  paragraphs. 

Figure  7-3  shows  the  schematic  of  the  flow  field  used  as  the  first  test 
example.  Initially  at  x  =  0  and  y  >  0  the  flow  is  inclined  at  15°  to  the  x 

axis.  At  x  =  0  and  y  =  0  the  flow  inclination  is  10°.  The  flow  has  a  Mach 

number  of  2  at  mesh  points  2-10  and  4  at  mesh  points  11-50  and  2.19  at  mesh 
point  1.  The  flow  situation  develops  as  follows.  There  is  initially  a 
discontinuity  in  total  pressure  between  mesh  points  10  and  11.  A  5°  expansion 
which  is  initially  between  mesh  points  1  and  2  propagates  into  the  flow  and 
interacts  with  the  contact  surface.  At  the  contact  surface  both  pressure  and 
flow  deflection  are  continuous  and  there  is  a  jump  in  stagnation  pressure. 

The  flow  is  divided  into  five  regions  by  the  expansion  waves  as  shown  in 
Figure  7-3.  In  these  regions  the  flow  properties  are  constant  and  the  exact 
solutions  are  given  in  the  table  on  the  fiqure.  Figure  7-4  shows  the 

numerical  computation  of  this  flow  using  KAYT1.  The  step  size  in  the  x 

direction  is  held  fixed  at  h/2  for  these  calculations.  The  values  at  the 

boundaries  y  =  0  (J  =  1)  and  y  =  49  (J  =  50)  are  also  held  fixed. 


Figures  (7-4)  (a)  and  (b)  show  the  smooth  development  of  pressure  and  flow 
deflection.  In  these  diagrams  it  is  hard  to  see  where  the  discontinuity  is. 
Figure  (7-4)  (c)  shows  the  total  pressure  profiles.  The  ramp  profile  joins 
the  constant  region  on  the  left  to  the  higher  but  constant  region  on  the 
right.  The  jump  in  total  pressure  takes  place  with  a  single  intermediate  mesh 
point.  At  two  stations  step  60  and  step  90  these  are  actually  two 
intermediate  mesh  points.  This  occurred  because  the  detection  scheme  for 
locating  the  total  pressure  discontinuity  was  very  primitive.  In  any  event 
the  algorithm  produces  self  mending  results.  Figures  (7-4)(d)  and  (e)  show 
the  profiles  of  Reimann  i  nvariants  9+  v  and  9—  »  respecti  vely.  The  code  does 
not  use  these  variables  in  the  computations.  These  graphs  serve  as  checks  on 
the  calculation  procedure.  The  jumps  that  occur  in  each  of  these  profiles  is 
a  result  of  the  jump  in  total  pressure  which  is  associated  with  a  jump  in  Mach 
number  and  hence  a  jump  in  v  the  Prandtl  Meyer  angle.  In  (7-4)(e)  the  initial 
expansion  fan  starts  at  x  =  0  between  mesh  points  1  and  2.  It  moves  to  the 
right  (along  9  + /j. characteristics )  and  spreads  out.  The  natural  spreading  of 
the  fan  is  augmented  in  the  numerical  calculation  because  the  wave  front 
(discontinuous  jump  in  derivative)  at  the  head  and  tail  of  the  fan  are  not 
tracked  by  the  program.  (There  are  no  general  algorithms  for  the  computations 
of  supersonic  flow  which  address  this  point.)  By  step  20  (x  =  10)  the 
expansion  has  intersected  with  the  contact  discontinuity  and  is  subsequently 
transmitted.  Only  a  wave  moving  to  the  right  appears  in  the  results  ford-*. 
This  is  in  accordance  with  classical  supersonic  flow  theory.  In  Fig.  (7-4)(d) 
the  left  moving  waves  become  evident.  Before  step  40  only  the  movement  of  the 
discontinuity  is  evident.  Then  at  step  40  the  reflected  expansion  wave  starts 
to  form.  The  values  of  9  +  v  and  6  -  »  in  regions  4  and  5  are  compared  with 


their  exact  values  in  the  table  below.  Is  it  clear  that  the  numerical 


calculation  is  very  accurate  in  the  computation  of  these  quantities. 


Comparison  of  Characteristic  Calculation  and  Finite  Difference  Calculation 


Region  f»^)Char  «>-»)FD 

4  .770  .771  -  .373  -  .373 

5  1.41  1.41  -1.01  1.01 


Figure  7-5  compares  the  computations  using  the  ramp  profile  with  a 
computation  that  does  not  specifically  “fit"  the  discontinuity.  Figure 
(7-5) (a)  compares  the  profiles  of  total  pressure.  The  profiles  are  smooth  but 
the  jump  in  total  pressure  is  spread  out  over  approximatley  fifteen  mesh 
points.  Figure  (7-5) (b)  and  (c)  show  the  profiles  of  the  Reimann  invariants. 
Again  the  results  are  smooth  and  have  the  correct  values  upstream  and 
downstream  of  the  discontinuity.  However,  it  is  clear  that  in  these  results 
the  discontinuity  zone  is  now  fifteen  mesh  points  wide.  This  represents  about 
one  third  of  the  mesh  points  used  in  the  calculation. 

The  most  series  errors  in  the  computed  results  using  the  ramp  function 
appears  in  figure  (7—4 ) (c ) .  The  exact  location  of  the  discontinuity  and  the 
center  of  the  ramp  function  have  drifted  apart.  This  is  a  result  of  the  wave 
front  spreading  noted  in  the  previous  discussion.  The  expansion  zones  spread 
out  ahead  of  the  exact  characteri sties  location  and  cause  the  value  of  the 
streamline  slope  to  change  in  advance  of  the  proper  positon.  In  order  to 


(a)  TOTAL  PRESSURE 

Fig.  7-5  Comparison  of  Numerical  Calculation  of  Test  Case  One 
With  and  Without  Ramp  Profile 
(Sheet  1  of  3) 


Fig.  7-5  Comparison  of  Numerical  Calculation  of  Test  Case  One 
With  and  Without  Ramp  Profile 
(Sheet  2  of  3) 


determine  if  this  is  the  correct  explanation  an  additional  calculation  was 
performed  that  did  not  have  a  wave  front  in  the  initial  profile.  Figure  7-6 
shows  the  characteristic  mesh  for  this  computation.  At  step  zero  there  is  an 
expansion  zone  which  is  propagating  toward  a  constant  pressure  boundary.  In 
the  finite  difference  program  this  constant  pressure  boundary  is  modelled  by  a 
jump  in  total  pressure.  The  expansion  zone  impinges  on  the  boundary  and 
reflects  as  compression  waves  which  move  back  into  the  flow.  On  this  figure 
is  the  computed  location  of  the  flow  boundary  by  the  finite  difference 
calculation.  The  agreement  is  virtually  perfect.  Figure  7-7  shows  pressure 
profiles  for  the  entire  computation.  Figure  7-7 ( a )  is  the  initial  pressure 
distribution  at  x  =  0.  Figure  ( 7-7 ) ( b)  is  at  x  =  10  where  the  expansion  has 
partially  reflected  from  the  boundary.  The  finite  difference  computation  and 
the  characteristic  calculation  are  in  virtual  agreement.  The  largest  errors 
seems  to  be  near  the  boundary  where  there  is  a  break  in  slope  in  pressure. 
Figures  7-7 (c)  and  (d)  shows  the  subsequent  development  of  the  flow  field.  At 
x  =  70  the  compression  wave  is  about  to  reflect  from  the  inner  boundary. 

There  is  a  spreading  of  the  wave  that  is  first  evident  at  x  =  40.  At  x  =  70 
the  finite  difference  results  first  preceded  then  lag  the  exact  calculation  by 
3-4  mesh  points.  The  precise  cause  of  this  dispersion  has  not  been  studied. 
Two  primary  sources  that  should  be  investigated  are  the  first  order  nature  of 
the  computational  algorithm  and  the  details  of  the  reflection  process  at  the 
boundary  ramp  function. 


Fig.  7-6  Test  Case  Two  ■  Calculation  of  Constant  Pressure  Boundary 


1.0 


Fig.  7-7  Comparison  of  Prassura  Profiles  for  Characteristic  Calculation 
and  Finite  Difference  Calculation  for  Test  Case  Two 
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The  use  of  the  ramp  function  in  three  dimensional  flow  requires  further 
work.  In  two  dimensional  flow  the  actual  position  of  the  discontinuity  was 
associated  with  the  intermediate  value  of  the  function.  In  essence  the  method 
tracks  the  discontinuity  on  a  subgrid  basis.  In  three  dimensions  this 
tracking  must  take  place  in  two  directions  (in  the  cross  plane).  A  method 
must  be  worked  out  to  accomplish  this  in  a  simple  manner.  The  ramp  function 
must  also  be  tested  out  in  flows  with  shock  waves.  The  shock  waves  have  a 
jump  in  all  flow  properties  and  not  just  a  single  one  as  was  employed  in  the 
test  calculations.  The  shock  propagation  problem  has  the  added  complication 
that  the  properties  behind  it  are  a  function  of  the  shock  (wave)  speed. 
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SECTION  8 


MACH  DISC  ANALYSIS  AND  NUMERICAL  COMPUTATION 

Shock  wave  systems  in  exhaust  plume  flow  fields  can  develop  in  a  manner 
which  produces  normal  shock  waves  or  Mach  discs.  In  axi symmetric  flows  the 
shock  wave  strength  intensification  process  leading  to  the  formation  of  the 
Mach  disc  is  understood  on  theoretical  grounds.  For  three  dimensional  plumes 
there  is  little  or  no  knowledge  of  the  processes  and  flow  properties  that 
result  in  normal  shock  segments  in  the  plume  (see  Section  3).  The  numerical 
prediction  of  the  location  of  the  Mach  disc  in  axisymmetric  flows  has  been 
accomplished  in  recent  years  (Refs.  8-1,  8-2).  These  studies  are  purely 
numerical  calculations.  In  this  section  an  analysis  is  derived  which 
explicitly  highlights  the  key  properties  of  the  exhaust  plume  that  influence 
the  Mach  disc  location.  Employing  the  results  of  this  analysis  a  new  numerical 
integration  scheme  is  derived  to  compute  Mach  disc  streamtubes.  This  new 
method  has  the  potential  for  overcoming  the  shortcomings  of  the  previous 
integration  procedures.  The  improved  understanding  of  the  axisymmetric  Mach 
disc  physics  should  serve  as  a  first  step  in  understanding  the  three 
dimensional  plume  situation. 

A  schematic  of  an  axisymmetric  plume  with  a  Mach  disc  is  shown  in  Fig. 
8-1.  At  all  points  in  the  flow  outside  of  the  Mach  disc  streamtube  the  plume 
gas  velocity  is  supersonic.  A  barrel  shock  divides  the  plume  into  a  "core" 
and  a  shock  layer  upstream  of  the  Mach  disc.  At  some  point  along  the  barrel 
shock  there  is  a  triple  point  which  is  the  confluence  of  the  barrel ,  reflected 
and  Mach  disc  shock  waves.  Downstream  of  the  triple  point  the  fl ow  i s  divided 
into  a  supersonic  flow  bounded  on  the  outside  by  the  reflected  shock  wave  and 
plume  boundary  and  on  the  inside  by  a  dividing  streamline  and  a  subsonic  core. 
The  subsonic  core  flow  has  passed  through  the  normal  shoe*  (Macb  disc'  and  >s 
centered  about  the  plume  axis. 


PLUME  BOUNDARY 


Fig.  8-1  Schematic  of  Axisymmetric  Underexpanded  Rocket  Exhaust  Plume 


According  to  Abbett  (Ref.  8-1)  the  triple  point  location  can  be  uniquely 
determined  by  finding  the  Mach  disc  position  which  leads  to  a  smooth  passage 
through  the  sonic  acceleration  of  the  subsonic  streamtube  (see  Fig.  8-1).  In 
both  Refs.  8-1  and  8-2  the  following  models  were  employed.  The  subsonic 
streamtube  was  described  by  a  one  dimensional  flow  model.  For  an  estimated 
position  of  the  triple  point  the  subsonic  flow  and  supersonic  flows  are 
computed  simultaneously  by  a  forward  integration  schemes.  Two  possible 
solutions  result:  a  supercritical  flow  has  the  subsonic  streamtube  reaching 
sonic  velocity  while  it  is  still  decreasing  in  area.  In  the  other  case,  a 
subcritical  flow  the  subsonic  streamtube  reaches  a  minimum  in  area  and  then 
increases  in  diameter  without  choking.  An  iterative  procedure  is  employed 
which  narrows  down  the  axial  position  of  the  triple  point  by  observing  these 
two  types  of  solutions.  This  is  a  classical  shooting  method  to  determine  a 
solution  for  a  system  of  equations  which  have  a  saddle  point  singularity 
downstream  (at  the  throat  of  the  subsonic  streamtube).  The  basic  disadvantage 
of  this  scheme  is  that  no  solution  for  the  subsonic  streamtube  emerges. 

Rather  a  bracket  on  the  triple  point  location  is  developed.  The  existence  of  a 
solution  within  this  bracket  is  not  certain. 

The  present  method  does  not  require  shooting  the  solution  from  the  triple 
point  location.  In  the  inviscid  approximation  the  location  of  the  sonic 
throat  can  be  calculated  without  a  subsonic  calculation.  The  subsonic 
streamtube  can  then  be  computed  by  an  upstream  integration  of  the  equations 
from  the  sonic  throat  making  use  of  certain  known  properties  of  supersonic 
flow  to  determine  the  outer  flow.  A  viscous  approximation  to  the  Mach  disc 
problem  is  presented  which  accounts  for  the  effects  of  mixing  between  the 
subsonic  streamtube  and  supersonic  flow  in  the  simplest  manner  possible.  This 


does  not  alter  the  basic  computational  scheme  but  adds  one  more  layer  to  the 
iteration  scheme. 

A  one  dimensional  model  for  the  subsonic  streamtube  is  retained.  The 
equations  governing  the  flow  passing  through  the  Mach  disc  are  (Ref.  8-3) 

M2-!  J_  dp_  +  dA_  1  +  (r-l)M2 
yM2  p  dx  A  dx  2 


1  OOsl 

p  dx 


(8-2) 


where  M  is  the  Mach  number,  p  the  static  pressure,  A  and  D  the  subsonic 
streamtube  cross  sectional  area  and  diameter,  p0  the  average  stagnation 
pressure,  and  'f'is  the  friction  factor.  Using  a  friction  factor  in  this 
manner  is  an  approximation  that  allows  for  the  simple  introduction  of  mixing 
effects  along  the  dividing  streamline.  The  flow  behavior  outside  the  dividing 
streamline  is  incorporated  in  the  function  f(x)  =  0+"  which  is  the  value  of 
the  Reimann  invariant  on  the  downward  running  characteristic  determined  from  a 
calculation  of  the  supersonic  flow.  The  pressure  on  the  outside  of  the  Mach 
disc  tube  flow  is  related  to  f  by  the  following  steps 

f  =  0'+  (8-3) 


(  )’  denotes  d/dx 

v  is  the  Prandlt  Meyer  function  and  is  only  a  function  of  external  Mach  number 


u  s  i  ng 

P  -y/y-l 

—  -  »  +4J  m2  > 

po2  2 

and  the  fact  that  p02  is  constant  along  the  bounding  streamline  equation 
(8-3)  becomes 


i  dp  =  ,f  ‘-e ' 
p  "dx"  h(m) 


(8-4) 


where  h(m)  =  p  dm  d>  =  -  1  ,, , 

U  ~T 


y-l  m2\  /  1  9 

^  ""T+€02  "  1+0 


2) 


where  0  =  and  €  =  (r-l)/(y+l) 


Equation  (8-4)  relates  the  pressure  gradient  in  the  external  supersonic  stream 
to  the  geometry  (  @'(x))  of  the  subsonic  streamtube.  The  solution  of  the  Mach 
disc  flow  requires  that  the  pressure  along  the  dividing  streamline  is  matched. 
Therefore,  combining  equations  (8-4)  and  (8-1)  the  basic  interaction  equation 
i  s 

i.j2-l  (f  '-9')  +  1  dA  =  1+CY-pM2  (4f)  (8-5) 

vm2  h(m)  A  dx  2  D 


For  a xi  symmetric  flow  A  =  tt y ^  where  y  is  the  local  height  of  the  subsonic 
streamtube.  The  governing  equations  become 


3-5 


“an* 


e‘  -f '  +  Be  =  a/  (i-M2) 

(3-6) 

-  ^3-  .  „ 

PQ  dx  ^ 

(8-7) 

where  B  =  ) 

(8-8) 

1-M2  y 

a  =  ?M2h(m)  [  l+(y-i)M2  ]  f/y 

(3-9) 

»i  ’  ■  yM2  .  (8-10) 

Equation  (3-6)  is  the  fundamental  interaction  equation  governing  the  pressure 
balance  between  the  supersonic  and  subsonic  flows  and  combined  with  equation 
(8-7)  (when^'^O)  must  be  solved  to  determine  the  subsonic  streamtube 
solution.  The  aoproximat ion  tan  Ofv  0  has  been  used  in  the  derivation.  The 
numerical  solution  of  equations  (8-6)  and  (8-7)  will  be  discussed  in  later 
paragraphs.  A  discussion  of  the  properties  of  equation  (8-6)  brings  out 
interesting  properties  of  the  Mach  disc  and  associated  plume. 

In  the  inviscid  case  it  is  possible  to  locate  the  sonic  point  in  the 
subsonic  flow  before  solving  equation  (8-6)  by  the  following  process,  for  a 
given  location  of  the  triple  point  along  the  barrel  shock  the  computation  of 
the  external  supersonic  flow  can  be  continued  employing  an  approximate  shape 
for  the  dividing  streamline  (see  Fig.  8-1).  From  this  solution  the  function 
f ( x )  =  8  +  v  along  the  dividing  streamline  can  be  determined.  Downstream  of 
the  Mach  disc  the  total  and  static  pressure  as  well  as  the  Mach  number  M| 
are  known  and  denoted  p^  and  p0j.  The  corresponding  sonic  pressure  p*  can 
be  computed  by  the  formula.  y 

p*  r  m  I  ~ 

?!  1  +  (V-DMj2 
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The  sonic  condition  requires  that  0=0.  Determine  the  pressure  on  the 
supersonic  side  of  the  dividing  streamline  at  every  point  for  0=0  from  the 


(known)  function  f(x) 

v  (x)  =  f  (x)  -  0  (x) 

"  (x)  =  f  (x) 

Use  the  isentropic  expressions  relating  »  to  m  and  m  to  p  to  compute  the 
pressure.  Figure  3-2  shows  a  schematic  plot  of  this  pressure  versus  x.  The 
sonic  point  or  throat  for  the  inner  flow  is  determined  at  the  point  where 
these  two  curves  cross.  (Figure  8-2  shows  two  such  crossings.  It  will  become 
clear  in  the  following  discussion  why  the  second  crossing  is  the  appropriate 
point).  The  point  denoted  x*  is  the  throat  because  it  has  the  two  required 
properties  6  =0  and  p  =  p*. 

Properties  of  the  solution  at  the  initial  point  and  the  throat  can  be 
deduced  from  the  interaction  equation.  For  inviscid  flow  the  interaction 
equation  is 

0'-f '+  B  0=  0  (3-11) 

The  coefficient  B  is  always  negative  for  subsonic  flow  (Note  that  h(m)  <  0). 

At  the  triple  point  the  pressure  is  matched  at  p^.  There  are  two 
possibilities  0  i  >  0  or  0  i  <  0.  For  0^  >0  0'-f ' =  -B0>  0  therefore 
0'>  f.  Thus  if  f'>  0  the  initial  dividing  line  will  be  unstable  because 
0'  >  0  and  the  slip  line  slope  increases  and  forms  an  upward  cusp.  To 
attain  a  concave  downward  curve  the  inital  value  of  f'  must  be  negative.  In 
terms  of  pressure  this  is  a  compression  wave.  In  the  case  of  0j  <  0  a 
stable  curve  can  only  be  achieved  for  f  <  0.  In  all  plume  calculations  that 
have  been  examined  the  initial  flow  behind  the  triple  point  is  compressive  or 
has  f'  >  0,  therefore,  in  these  cases  it  is  concluded  that  the  initial  slope 
of  the  dividing  streamline  must  be  positive. 
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f  <  0 
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Fig,  8-2  Locating  the  Sonic  Throat  by  Matching  Supersonic 
Side  Pressure  to  Sonic  Pressure 
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At  the  sonic  throat  the  flow  is  accelerating  and  so  dp/dx  <  0.  The 
pressure  gradient  has  the  same  sign  as  9'- f'  since  h(m)  <  0  (see  equation 
(8-4)).  Thus  9'-f‘  <  0  or  f1  >0'.  The  sonic  throat  has  the  following  two 
conditions  9  -  0  and  B'  >  0.  Therefore,  the  throat  can  only  occur  when 
f'  >  0  or  an  expansion  portion  of  the  supersonic  flow.  This  precludes  the 
sonic  throat  from  occurring  at  the  first  place  the  curves  cross  in  Fig.  8-2 
where  f'<  0. 

The  following  general  picture  (Fig.  3-3)  emerges  for  the  Mach  disc 
streamtube.  The  solution  curve  in  Fig.  8-3  is  the  pressure  on  the  dividing 
streamline  or  in  the  subsonic  streamtube.  Also  shown  is  the  curve  of  pressure 
in  the  supersonic  stream  for  9  =  0  on  the  dividing  line  (denoted  curve  A)  and 
the  values  of  total  pressure  p0^  and  sonic  pressure  p*  in  the  subsonic 
streamtube.  When  the  solution  curve  is  above  (below)  curve  A  9  >  0  {9  <  0). 
The  solution  starts  downstream  of  the  triple  point  with  pressure  p^.  For 
all  flows  where  there  is  a  compression  following  the  triple  point  in  the 
supersonic  flow  the  initial  angle  must  be  positive  and  decreases  with  x. 
Downstream  of  the  triple  point  the  subsonic  streamtube  undergoes  an  increase 
in  pressure  until  the  value  of  9  reaches  zero.  Downstream  of  this  point  in 
Zone  (2)  f'  is  still  negative  signifying  a  compressive  outer  wave  pattern  but 
O'  <  f'  so  that  the  pressure  drops  and  the  angle  decreases.  In  Zone  3  the 
outer  flow  is  expanding  and  the  function  f1  >  0.  The  value  of  <9  increases 
smooth1y  (  9'  >0)  and  reaches  zero  precisely  at  the  point  where  sonic 
conditions  are  met  and  the  inner  streamtube  is  choked. 
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Fig.  8-3  Geometry  of  Mach  Disc  Stream  Tube  and  Corresponding 
Solution  of  Interaction  Equation 
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It  is  clear  that  the  function  f  plays  a  decisive  role  in  determining 
the  subsonic  streamtube  solution  and  hence  Mach  disc  position.  The  value  of 
the  function  f  is  related  in  a  complex  way  to  the  entire  plume  flow  solution 
upstream  of  the  triple  point.  It  is  possible,  however,  to  highlight  some 
interesting  properties  which  determine  f.  The  value  of  f,  the  Reimann 
invariant  on  the  downward  running  characteristic,  would  remain  constant  and 
equal  to  its  value  on  the  plume  boundary  if  the  flow  were  two  dimensional  and 
isentropic.  Because  the  plume  is  axi symmetric  and  noni sentropic  f  is  not 
constant  along  the  waves  but  it  is  expected  that  'whatever  variations  are  due 
to  these  effects  are  approximately  the  same  for  all  waves. 

Figure  8-4  is  a  schematic  of  two  plumes  at  widely  different  free  stream 
conditions.  The  quiescient  plume  is  shown  for  a  slightly  underexpanded  nozzle 
exit  conditions.  The  Mach  number  in  the  plume  shock  layer  (between  interface 
and  barrel  shock)  is  low,  therefore,  the  downward  running  waves  are  at  a  steep 
angle  and  reach  the  Mach  disc  streamtube  in  a  short  axial  distance.  The 
pressure  on  the  plume  boundary  is  constant  so  that  the  variation  in  the 
function  f  comes  from  the  curving  plume  boundary  (Constant  pressure  means  that  v 
is  constant).  The  extent  of  the  plume  interface  which  directly  influences  the 
Mach  disc  flow  is  rather  short.  The  highly  underexpanded  plume  (Fig.  8-4b) 
has  quite  a  different  picture.  The  gas  in  the  plume  shock  layer  is  at  very 
high  Mach  number  so  that  the  portion  of  the  plume  boundary  where  f  originates 
is  quite  extended.  In  addition  the  plume  is  in  a  hypersonic  free  stream  which 
creates  a  large  pressure  gradient  on  the  interface.  Therefore,  the  details  of 
the  pressure  and  deflection  along  a  large  segment  of  the  plume  boundary  are 
influential  in  determining  both  9  and  v  along  the  interface  and  hence  on  the 
value  of  f  at  the  Mach  disc  interface. 
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PRIMARY  ZONE  OF  INFLUENCE 


(a)  QUIESCENT  BACKGROUND 


lb)  HIGHLY  UNDEREXPANDED  PLUME  IN  A  HYPERSONIC  STREAM 

2193-OS8D 

Fig.  8-4  The  Boundary  Zone  of  Influence  on  the  Mach  Oise  Flow 
for  Two  Types  of  Flow  ' 
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The  effect  of  viscous  mixing  along  the  Mach  disc  streamline  introduces 
additional  complexity  and  new  features  in  the  solution.  Within  the  scope  of 


one  dimensional  analysis  it  is  possible  to  introduce  mixing  effects  in  a 
global  sense.  It  is  realistic  to  employ  this  approximation  at  present  because 
the  actual  detailed  two  dimensional  flow  field  is  a  complex  transonic  flow 
with  turbulent  mixing  and  transverse  pressure  gradients  which  would  be 
difficult  to  accurately  model.  In  the  model  adopted  here  the  influence  of  the 
mixing  is  to  introduce  a  mechanism  that  increase  the  total  pressure  in  the 
subsonic  flow.  Mixing  also  changes  the  flow  inclination  that  the  supersonic 
flow  sees  along  the  dividing  streamline  through  a  displacement  effect.  This 
effect  is  not  incorporated  in  the  present  model.  The  sonic  throat  condition  is 
no  longer  9  =  0  but  it  has  a  small  negative  value.  The  increasing  value  of 
p0  with  distance  results  in  increasing  flow  Mach  number  even  at  constant 
static  pressure.  Equations  (8-6 ) - ( 8-10 )  show  the  relationship  of  the  mixing 
terms  in  the  interaction  equation  (3-6)  and  the  total  pressure  equation  (8-7). 
The  Mach  disc  flow  has  lower  velocity  and  lower  total  pressure  than  the  outer 
supersonic  flow  so  that  f  is  negative  and  cr  i  is  positive.  A  solution 
curve  with  mixing  is  sketched  in  Figure  8-5.  The  main  difference  to  note  in 
this  solution  compared  to  the  inviscid  solution  (Fig.  8-3  )  is  that  p0  and 
p*  are  increasing  with  distance.  The  sonic  throat  is  again  located  at  the 
point  where  the  pressure  (ps)  evaluated  by  setting  9=9*  (using  the 
function  f)  is  equal  to  p*.  In  this  case,  however,  the  value  of  9  *  is  not 
known  a  priori  and  must  be  determined  as  part  of  the  solution. 
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The  numerical  integration  procedure  for  the  so’ut’on  of  tne  ’nteract’on 
equation  is  based  on  integrating  the  equation  (3-61  oackward  starting  at  tne 
sonic  throat.  Several  integration  schemes  have  been  employed  with  virtual 'y 
equal  success.  The  numerical  integration  procedure  for  the  'nviscid  case  's 
outlined  first.  The  extension  to  the  viscous  case  follows  naturally  and  is 
detailed  at  the  end  of  this  chapter.  The  starting  solutions  at  the  throat  are 
derived  and  their  numerical  implementation  is  discussed,  'he  solution 
proceeds  using  the  following  sequence  of  steps: 

(In  the  follownq  discussion  superscript  (i)  denotes  iteration  number  and 
yO)  denotes  the  subsonic  streamtube  height  computed  as  a  result  of  the 
estimated  value  by  yO).) 

1)  Compute  y*  and  estimate  M  ( * ) ,  yU) 

2)  Determine  the  location  of  the  throat  x*  and  the  flow  gradients  there. 

3)  Solve  the  interaction  equation  for  9  (x)  by  integrating  backwards 
from  the  point  x*,  y*  where  0=0 

4)  Determine  yO)  by  integrating  the  above  solution  for  9  . 

5)  Combine  the  estimated  values  of  yO)  and  the  computed  values  of  Y^) 
to  achieve  the  next  estimate  for  the  Mach  disc  streamtube  y^+^) 

6)  Repeat  the  steps  (3)  -  (5)  until  the  solution  converges. 

The  details  for  each  step  are  outlined  in  the  following  paragraphs. 

Step  2  The  solution  of  the  interaction  equations  has  a  classical  saddle  point 
singularity  at  the  sonic  throat.  This  can  be  verified  by  assuming  that  the 
Mach  number  in  the  neighborhood  of  the  throat  is  given  by 
M^  =  1  +  a  (x  -  x*)  and  investigating  the  nature  of  the  interaction  equation 
near  x  =  x*. 
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lonsiaenng  the  relationsmp  between  press  ure  ana  ^acn  number  ‘or  tne  Macn 


Use  flow 


_1 _ dM“  =  -2  _1  d£ 

EM2  dx  p  lx 

E  =  1  +il!’*2 


O  i  up  _  ■ n  ■  -  •  .  -  . 

Substitute  =  1  *  a  (x  -  x’}, - r  \9  -t  ;  to  i  m ;  and  ^  ^  o'x  -x*'  to 

p  dx 

get  a  relationship  between  the  Mach  number  gradient  a  and  tne  dividing 

streamline  angle  gradient  b. 
y 

b  •=  -  — -  h*  a  +  f  * 

y  ri 

In  order  to  be  self  consistant  the  subsonic  flow  must  also  satisfy  the  one 


dimensional  Mach  number  area  relationship 


1-M2  dM2 
EM2  dx 


2  dA  -46 
A  dx  y 


Again  using  the  expansions  for  M  and  9  near  the  throat  a  second  relationship 
between  a  and  b  can  be  deri ved 


Combining  these  two  relationships  results  in  a  single  quadratic  equation 


ay*  +  a  (2yh*)  =  2  ( y  +  1)  f  ’* 

Step  3  The  interaction  equation  (8-6)  is  integrated  from  one  node  point  to 
the  next.  The  value  of  the  function  B  is  assumed  known  from  a  previous 
iteration  step.  If  i  denotes  the  iteration  step  and  x  1  is  the  value  of  x 


for  mesh  poi  nt  j 


e‘-f  '  +  Bfl  =  o 


e'+1  -  eul  -  <Vi  '  V  *  f3' 
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The  evaluation  of  the  integral  in  equation  (8-12)  was  achieved  by  either  of 
two  methods  both  of  which  are  given  below.  In  method  1  the  trapezoidal  rule 


is  employed  resulting  in  ...  ■ 

(1  +B7  Ax/2)  -  (f.  .  -  f.) 

„i+:  3  ^  M  J+1  J 

(1-B.7Ax/2) 

j 

In  the  second  method  it  is  assumed  that  the  function  B  has  a  singular  behavior 
as  fol lows 

8*  a  ,  , 

3  ( x )  =  ~ —  +  3  ^  x ) 


where  8*  is  a  constant  which  characterizes  the  singularity  at  x*.  This 
results  in  an  integration  formula 

<;  !'•  *»/2]  -  :V!-M 

1  -  ,3  3  1  -  Q.1  Ax/2 


where  a  - 


x  ^ ;  /  v  x  -  ,  -  x;!  v  h 

J  J  *  u 


(  X  •  - 

J 

KJ- 

^  * 

j  1  <■ n  'X  *  Xj 

)/(  x 

fprmu'la  is  used 

to  i 

t  1  -  !**  -  x  ,•  );'(*. 


backwards  the  point  x  =  <*,  y  =  y*  usinq  the  solution  for  9 

-r-  =  tan  t?  «  9 


; .  .  -  x  . ) 

J+1  j 


lenot. 1  ng  the  va'ue  by  V 


f.  1  =  1 


'  -*■  1  1  ♦  ' 
^  '  r  +  h  . 

:  ti  j 


Step  b  Based  on  tne  new  value  of  d,'*-  at  each  mesh  point  an  updated 

/a I  ue  ‘or  the  ‘^ach  number  on  the  supersonic  side  of  the  l^vdinq  streamline 

i-l 

tan  oe  -.omputed  by  invert!  ho  •  for  T;'*'  js'na 


Using  this  Mach  number  and  the  total  pressure  of  the  outer  (supersonic)  stream 
an  updated  value  of  static  pressure  can  be  computed.  Combining  this  pressure 
and  the  total  pressure  of  the  subsonic  streamtube  at  each  point  with  the  Mach 
number  pressure  relationship  yields 

sr={tT[(v^-f 

Based  on  the  computed  value  Yji  +  1  and  the  known  value  of  y*  for  the 

subsonic  streamtube  another  estimate  of  the  subsonic  streamtube  Mach  number  is 

computed 

’Mji+1  =  F  (Yji+l/y*) 

where  the  function  F  is  used  symbolically  to  denote  the  inverse  of  the  Mach 
number  -  A/A*  relationship.  (This  equation  is  actually  solved  by  a  simple 
iteration.)  Finally  a  new  value  of  M  is  determined  by  combining  the  estimates 
M  and'M'  using  the  underrelaxation  formula 

Mj i *1  =  .25  +  1  +  -75  (Mj  +  *j)/? 

The  solution  for  the  viscous  case  proceeds  in  exactly  the  same  manner  as 
outlined  above  for  the  inviscid  streamtube  case  with  the  following 
modifications.  The  initial  location  of  the  sonic  point  cannot  be  made  because 
the  sonic  throat  angle  9 *  is  given  by  the  formula 

9*  --  yf*/2 

Since  "f1*  depends  on  the  solution  an  initial  estimate  of  the  viscous  solution 
is  necessary  to  determine  x*.  The  equation  for  stagnation  pressure  p0  of 
the  inner  streamtube  is  achieved  by  straightforward  integration,  again  based 
on  the  current  iterate  value  of  a i  (x).  In  the  inversion  to  determine 

y\ 

Mji+1  the  value  of  y*  must  be  computed  for  each  station  using  the 
relationship  p0A*  =  constant  for  tne  inner  flow  . 
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Two  test  cases  were  employed  to  determine  if  the  theory  and  computer 
code  were  operating  satisfactorily.  The  first  test  case  consisted  of  a  known 
exact  solution  (inviscid)  which  the  program  was  required  to  reproduce.  The 
second  test  case  was  an  actual  Mach  disc  flow  with  substantial  mixing  effects. 
In  the  first  case  the  computer  code  and  theory  reproduced  the  exact  solution 
to  a  very  close  tolerance.  In  the  second  case,  the  code  was  used  to  analyze 
the  flow  pattern  and  achieved  very  good  agreement  with  all  reported  measured 
quantities.  However,  the  experimental  results  reported  in  Reference  (8-4) 
were  not  complete  enough  to  allow  prediction  of  the  Mach  disc  location. 

The  first  test  case  was  devised  by  prescribing  a  subsonic  streamtube  given 
by  the  equation 

y  =  1.36  +  .36  cos  (.59169x  -  2.1664)  (8-13) 

which  is  shown  in  Figure  8-6.  The  streamtube  is  assumed  to  choked  at  the 
throat  where  y  =  1.0  the  Mach  number  is  unity.  Using  isentropic  flow  tables 
the  Mach  number  distribution  and  pressure  distribution  in  the  streamtube  were 
computed  .  The  compatible  outer  supersonic  flow  is  constructed  as  follows.  A 
total  pressure  ratio  (outer  supersonic  stream  to  inner  subsonic  stream) 

P02/P0I  =  l4- 403  was  chosen.  The  outer  static  pressure  is  set  equal  to 
the  inner  subsonic  streamtube  pressure.  The  pressure  and  total  pressure  is 
known  at  each  point  on  the  slip  line  so  that  the  Mach  number  can  be  computed 
( m( x ) ) .  The  local  slope  of  the  streamtube 

Q  =  y'  =  -.21301  sin  (.59169x  -  2.1664) 

is  computed  at  each  point  and  thus  the  function  f(x)  =  9  +  »  can  be  computed 
on  the  supersonic  side  of  the  slip  line.  The  Prandlt  Meyer  function  ^  is  only 
a  function  of  m  and  is  computed  at  each  point.  This  construction  has  defined 
a  subsonic  streamtube  with  a  consistant  value  of  9  +  v  resulting  in  matched 
pressure  on  the  dividing  streamline. 
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Based  on  the  following  input  f(x),  p02/Pol  =  14.403,  m(x  =  0)  =  .5, 
y  (x  =  0)  =  1.158  the  computer  code  is  used  to  predict  the  slip  line  height 
y(x).  If  the  program  is  working  properly  it  should  reproduce  the  streamtube 
given  by  equation  (8-13).  The  computed  solution  can  then  be  compared  with  the 
original  equation  to  determine  the  accuracy  and  efficiency  of  the  computer 
code.  Tables  1  and  2  summarize  these  comparisons.  Table  1  compares  the  exact 
original  geometry  and  the  computed  numerical  results.  Columns  three  and  four 
are  the  errors  in  the  streamtube  height  and  slope  respectively.  This 
calculation  employed  forty  iterations  and  used  method  2  for  the  integration 
scheme.  The  errors  are  at  most  2  x  10“3  which  is  entirely  acceptable  in 
this  case.  (The  function  f(x)  was  only  input  to  three  decimal  places.)  For 
this  calculation  the  pressure  on  the  subsonic  and  supersonic  side  of  the  slip 
line  are  compared  in  Table  2.  Column  one  is  the  accumulated  difference  (sum 
of  the  absolute  values)  in  pressure  difference  along  the  slip  surface.  The 
total  accumulated  error  is  .0073  for  the  entire  37  points  of  the  calculation. 
Thus  on  the  average  the  difference  in  pressure  at  each  mesh  point  is 
approximately  .0002.  This  is  well  within  the  expected  calculation  accuracy. 

The  convergence  of  the  solution  is  shown  in  Figure  8-7,  where  the  Mach 
number  of  the  subsonic  streamtube  is  plotted  for  the  iterations  1,2, 3, 5  and 
10.  Iteration  1  is  the  starting  solution  which  in  this  case  was  chosen  to  be 
M  =  .25  everywhere  but  in  the  vicinity  of  the  throat  where  the  local  throat 
solution  was  employed  to  get  a  good  first  approximation.  The  relaxation 
procedure  showed  very  good  properties  in  this  case.  The  successive  estimates 
for  Mach  number  proceeds  smoothly  (but  not  monotonically)  toward  the  final 
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Table  1  Comparison  of  Numerical  Calculation  and  Exact  Solution 


J 

X 

y~yexact 

y ‘ -  y ' exact 

^exact 

y' exact 

xlO2 

xlO2 

1 

.00 

-.0949 

-.0176 

1.158037 

.17633 

2 

.25 

-.0985 

.0564 

1.204165 

.19202 

3 

.50 

-.0648 

.2847 

1.253695 

.20351 

4 

.75 

-.0493 

-.0842 

1.305549 

.21056 

5 

1.00 

-.0706 

-.0096 

1.358591 

.21301 

6 

1.25 

-.1015 

-.1594 

1.411664 

.21080 

7 

1.50 

-.1082 

.1811 

1.463609 

.20400 

8 

1.75 

-.1058 

-.0888 

1.513290 

.19273 

9 

2.00 

-.1107 

.1167 

1.559624 

.17726 

10 

2.25 

-.1322 

-.2269 

1.601598 

.15792 

11 

2.50 

-. 1741 

-.  0551 

1.638295 

.13512 

12 

2.75 

-.1880 

-.0112 

1.668913 

.10938 

13 

3.00 

-.1873 

.0519 

1.692785 

.08125 

14 

3.25 

-.1744 

.0754 

1.709388 

.05134 

15 

3.50 

-.1546 

.0955 

1.718360 

.02031 

16 

3.75 

-.1373 

.0452 

1.719505 

-.01116 

17 

4.00 

-.1234 

.0559 

1.712798 

-.04239 

18 

4.25 

-.1065 

.0588 

1.698386 

-.07270 

19 

4.50 

-.0968 

-.0139 

1.676582 

-.10141 

20 

4.75 

-.0969 

-.0273 

1.647865 

-.12791 

21 

5.00 

-.1112 

-.1430 

1.612859 

-.15162 

22 

5.25 

-.0992 

.1835 

1.572332 

-.17201 

23 

5.50 

-.0679 

.0013 

1.527166 

-.18865 

24 

5.75 

-.0381 

.1653 

1.478350 

-.20117 

25 

6.00 

-.0160 

-.0631 

1.426948 

-.20929 

26 

6.25 

-.0036 

.0819 

1.374085 

-.21285 

27 

6.50 

.0274 

.0859 

1.320913 

-.21175 

28 

6.75 

.0502 

.0201 

1.268596 

-.20603 

29 

7.00 

.0459 

-. 1282 

1.218274 

-.19581 

30 

7.25 

.0176 

-.1674 

1.171049 

-.18131 

31 

7.50 

.0186 

.1129 

1.127950 

-.16285 

32 

7.75 

.0556 

.  1274 

1.089919 

-.14084 

33 

8.00 

.0614 

-. 1276 

1.057787 

-.  11575 

34 

8.25 

.0619 

.0942 

1.032257 

-.08813 

35 

8.50 

.069? 

-.0625 

1.013884 

-.05859 

36 

8.75 

.0565 

-.0553 

1.003071 

-.02776 

37 

9.00 

.0397 

-.0854 

1.0C0053 

-.00367 
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Table  2  Comparison  of  Pressures  on  Either  Side  of  the  Slip  Surface 
J  ^  AP  x  102  Psubsonic  Psupersonic 

i  =  i 


1 

.0145 

.842616 

.842761 

2 

.0324 

.870053 

.870231 

3 

.0522 

.892628 

.892827 

4 

.0748 

.910797 

.911023 

5 

.0983 

.925232 

.925472 

6 

.1223 

.936721 

.936962 

7 

.1462 

.945849 

.946088 

8 

.1700 

.953059 

.953297 

9 

.1932 

.958691 

.958923 

10 

.2162 

.963049 

.963279 

11 

.2388 

.966365 

.966590 

12 

.2610 

.968865 

.969087 

13 

.2829 

.970659 

.970878 

14 

.3045 

.971839 

.972055 

15 

.3259 

.972460 

.972674 

16 

.3472 

.972546 

.972760 

17 

.3686 

.972106 

.972320 

18 

.3903 

.971118 

.971335 

19 

.4124 

.969525 

.969746 

20 

.4347 

.967241 

.967464 

21 

.4573 

.964141 

.964367 

22 

.4808 

.960097 

.960331 

23 

.5046 

.954893 

.955131 

24 

.5286 

. 948238 

.948479 

25 

.5526 

.939761 

.940001 

26 

.5764 

.929029 

.929268 

27 

.5993 

.915578 

.915807 

28 

.6206 

.898747 

. 898960 

29 

.6393 

.877762 

.877949 

30 

.6545 

.851860 

.852012 

31 

.6673 

.820728 

.820857 

32 

.6765 

.784111 

.784203 

33 

.6777 

.741203 

.741215 

34 

.6791 

.692629 

.692642 

35 

.6874 

.639292 

.639209 

36 

.7033 

.581848 

.581689 

37 

.7279 

.522517 

.522270 

_ I _ I _ I _ I _ L 

2  3  4  5  6 

Fig,  8-7  Convergence  of  Mach  Numbers  for  Inviscid  Test 
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result  which  are  shown  in  Figure  8-3.  In  this  particular  example  the  function 
f  actually  has  a  slight  increase  near  x  =  0  which  is  reflected  in  the  initial 
increase  in  slope  of  the  slip  line  (  0‘(O)  >  0).  Mach  number  decreases  from 
its  initial  value  of  .5  to  .2  at  x  =  4  with  the  associated  increase  in 
pressure  to  near  stagnation  pressure.  In  this  region  of  very  low  Mach  number 
the  subsonic  streamtube  behaves  very  much  like  a  constant  pressure  boundary. 

In  mathematical  terms  the  function  3  in  the  interaction  equation  is  very  small 
and  to  a  good  approximation  the  interaction  equation  reduces  to0-  f'  =  0  with 
the  solution 

9  =  f  -  f(  8=  0) 

This  property  is  displayed  graphically  in  Figure  8-9  where  9  and  the 
difference  f  -  f(0=O)  are  plotted.  In  this  type  of  flow  the  central  role  of 
the  function  f  is  very  clear.  Only  near  the  sonic  throat  where  3  is  order 
unity  does  the  solution  for  9  differ  markedly  from  f  -  f  (  9  -  0) 

The  second  test  case  employed  the  flow  field  in  a  supersonic  diffuser  as 
reported  in  Ref.  8-4.  The  geometry  of  the  flow  is  shown  in  Figure  8-10.  An 
incident  shock  approaches  the  diffuser  axis.  At  z  =  0  a  Mach  disc  is  formed  at 
the  intersection  of  the  reflected  and  incident  shock  waves.  There  is  a  mixing 
layer  between  the  supersonic  stream  and  the  flow  which  passed  through  the  Mach 
disc  (normal  shock).  In  the  mixing  layer  the  flow  Mach  number  smoothly  passes 
from  subsonic  to  supersonic.  Therefore,  the  sonic  line  appears  to  eminate 
from  the  triple  point  and  move  almost  horizontally  at  first  before  curving 
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Fig.  8-10  Flow  Field  for  Test  Case  2  •  Viscous  Mach  Disc 
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downward  and  intersecting  the  axis.  Reference  8-4  reports  the  profiles  of 
Mach  number,  total  pressure  and  impact  pressure  at  six  stations  downstream  of 
the  shock  triple  point.  Only  the  first  three  stations  are  relevant  to  this 
study  because  the  inner  streamtube  is  supersonic  downstream  of  these.  It  is 
important  to  note  that  the  exact  location  and  size  of  the  Mach  disc  were  not 
measured  in  the  experiment  and  that  the  closest  station  to  the  triple  point  is 
at  z  =  .5  inches.  It  is  pointed  at  in  reference  8-4  that  no  attempt  was  made 
to  probe  the  exact  for  orientation  of  the  shock  configuration  in  the  vicinity 
of  the  triple  point  where  the  shock  may  be  curved. 

This  flow  field  was  studied  using  the  present  analysis  as  follows.  Since 
the  details  of  the  upstream  flow  field  are  unknown  not  enough  information  is 
presented  to  attempt  an  a  priori  computation  of  the  Mach  disc  location  and 
height.  It  is  possible,  however,  to  determine  if  the  computer  code  gives 
results  consistant  with  various  aspects  of  the  experiment.  In  order  to  do  this 
the  function  f(x)  =  9  +  v  was  computed  using  the  experimentally  observed 
values  for  the  Mach  number  on  the  supersonic  side  of  the  slip  line  and  the 
reported  slip  line  geometry.  From  the  published  results  the  value  of  the 
streamline  slope  was  taken  to  be  approximately  zero.  This  was  all  possible 
downstream  of  the  first  measured  station  which  was  at  z  =  .5  inches.  The 
properties  at  the  triple  point  and  its  precise  location  were  not  presented  in 
the  reference.  The  triple  point  solution  was  constructed  subject  to  the 
constraint  that  the  pressure  downstream  was  given  by  the  normal  shock  pressure 
(M  =  4.60  upstream).  This  renders  the  Mach  disc  calculation  a  function  of 
only  a  single  parameter  -  the  deflection  across  the  incident  shock,  for 
example.  It  was  determined  that  the  Mach  number  M4  (see  Fig.  8-10) 
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downstream  of  the  triple  point  (supersonic)  was  only  a  weak  function  of  the 
initial  deflection  angle  (  6  )  and  equal  to  1.89  which  is  entirely  consistent 

with  the  experimental  results.  The  initial  flow  angle  downstream  of  the 
triple  point  (  6  .  see  Fig.  8-10)  is  a  strong  function  of  5^.  Since  the 
flow  field  is  expanding  downstream  of  the  triple  point  the  initial  angle  must 
be  negative.  The  value  #i  was  prescribed  to  be  -3.5°  and  the  solutions  were 
not  found  to  be  inconsistant  with  this  value. 

With  function  f  constructed  above,  the  value  of  P02/P0I  based  on  the 
triple  point  solution  and  an  estimated  initial  height  of  the  Mach  disc  a 
series  of  solutions  was  computed.  Only  one  series  of  computations  are 
detailed  here.  Three  values  of  the  friction  factor  constant  were  employed 
-.01,  -.005  and  0  (the  inviscid  case).  Figure  8-11  shows  the  computed  slip 
line  location  with  two  experimentally  measured  Mach  number  profiles  together 
with  the  result  of  the  computations  for  the  one  dimensional  results.  It  is 
important  to  note  that  the  slips  line  location  quoted  in  reference  8-4, 
denoted  by  S  are  far  outside  of  an  estimate  based  on  the  Mach  number  profile. 
It  is  believed  that  this  occurred  because  the  authors  extrapolated  the  stream¬ 
lines  back  to  a  triple  point  location  they  estimated  based  on  a  straight 
transmitted  shock  from  z  =  .5  back  to  z  =  0.  In  the  present  calculation  the 
initial  height  of  the  Mach  disc  was  chosen  so  that  the  slip  line  passed 
approximately  through  the  center  of  the  mixing  layer  (based  on  Mach  number 
profile)  at  the  z  =  .5  station.  At  both  stations  the  one  dimensional  values 
of  Mach  number  are  reasonable  averages  of  the  measured  profiles.  Mote  that 
the  one  dimensional  values  have  gone  supersonic  at  the  second  station. 


Fig.  8-11  Viscous  Test  Case  ■  Computed  Slip  Line  and  Mach  Numbers 


Figure  8-12  presents  the  computed  slip  line  location  for  the  three  cases 


and  the  Mach  number  calculations  along  with  the  experimental  centerline 
values.  The  effect  of  viscosity  is  to  move  the  throat  upstream  closer  to  the 
Mach  disc.  This  is  expected  because  mixing  in  these  cases  increase  the 
average  total  pressure  of  the  subsonic  stream  which  increases  the  Mach  numbers 
(by  decreasing  the  local  value  of  p/p0).  The  experimental  sonic  line  is 
shown  in  the  figure.  It  is  not  possible  from  these  results  to  conclude  which 
value  of  friction  factor  is  most  appropriate.  The  measured  centerline  Mach 
numbers  are  below  the  calculated  values  in  all  crses  which  is  to  be  expected 
based  on  the  upward  curvature  of  the  dividing  streamline  . 

Figure  8-13  shows  the  total  pressure  distribution  for  the  three  cases.  The 
total  pressure  increment  at  the  sonic  point  is  13  %  and  8  "  for  the  cases  with 
friction  constants  k  =  -.005  and  -.01  respectively.  In  this  case  the  effect 
of  mixing  does  not  have  a  dramatic  effect  on  the  geometry  of  the  flow  field 
(Fig.  8-12).  Friction  moves  the  sonic  point  a  noticeable  amount,  however,  by 
far  the  largest  force  on  the  flow  is  the  static  pressure  gradient.  These 
conclusions  cannot  be  carried  over  to  the  Mach  disc  in  a  plume  because  in  that 
case  the  initial  pressure  gradient  is  compressive  behind  the  Mach  disc  (f(x) 
decreases).  Therefore,  the  static  pressure  gradient  tends  to  decelerate  the 
flow  in  opposition  to  the  total  pressure  gradient  which  is  driving  Mach  number 
in  the  opposite  direction. 
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Fig.  8-12  Comparison  of  Computed  Results  for  Various  Friction  Factors 
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NOTES  FOR  VISCOUS  SOLUTIONS: 

Throat  solution  -  the  Mach  number  gradient  a  at  the  throat  is  g 
solving  the  equation 


a2(y*)  +  a(2rh*  -  (y+i)yf*)  =  2(y+i)f'*  -y(y+i)£‘*  +  2y(y+i)f* h*/y» 


The  slope  gradient  b  is  then 

b  =  f‘*  -yh*a/(y+l)  +yh*£*/y* 


and  the  throat  angle  is 

d*  =yf* 

The  friction  factor  f  is  related  to  the  inner  and  outer  Mach  numbers 


v i  + 


_  -  M 

- 1  J  i  + 
2 


2 


where  f  =  2r/pv^  and  k  and  Traf  are  parameters  in  the  eddy  viscosity 
for  tne  shear  layer. 
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SECTION  9 


CONCLUSIONS 

The  optical  signature  and  thermochemical  properties  of  rocket  exhaust 
plumes  are  sensitive  functions  of  the  many  parameters  determining  the  plume 
flow  field.  The  influence  of  the  many  inputs  bearing  on  the  plume  signature 
is  basically  through  the  temperature  and  species  fields.  The  plume  is  a  hot 
mass  of  gas  composed  of  a  variety  of  optically  active  species.  The  precise 
determination  of  the  optical  signature  is  based  on  accurate  fluid  mechanical 
predictions  leading  to  spatially  resolved  temperature,  pressure  and  species 
concentration  fields.  Both  chemical  kinetic  and  radiative  transport  processes 
are  driven  by  source  or  rate  terms  which  are  of  the  Arrenhius  type. 

Therefore,  the  IR  signature  is  very  sensitive  to  temperature  and  this 
sensitivity  increases  as  the  level  of  plume  temperature  decreases.  A  model 
plume  has  been  analyzed  based  on  mathematical  analysis  taking  advantage  of 
this  mathematical  property.  The  direct  quantitative  relationship  between 
temperature  field  and  species  field  and  local  station  radiation  has  been 
shown. 

In  the  general  case  of  multinozzle  rocket  plume  flow  fields  the  inviscid 
pattern  is  a  complex  three  dimmensional  flow  containing  several  shock  wave 
surfaces.  The  shock  waves  produce  both  near  and  far  field  temperature 
increases  and  so  are  central  to  (chemical  activity  and)  optical  signature 
predictions.  While  the  shock  wave  structure  is  not  the  only  fluid  dynamic 
process  involved  in  determining  the  temperature  distribution  it  is  involved  in 
inviscid  flow  field  calculation  which  is  the  primary  skeleton  on  which  is 
built  the  total  flow  field  picture.  Accurate  prediction  of  the  shock  wave 
structure  can  be  achieved  only  through  detailed  calculations  which  track  the 
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shock  surfaces.  Methods  which  "capture"  the  shock  waves  would  require  machine 
memories  and  computer  times  which  are  far  beyond  what  is  presently  available. 

A  floating  three  dimensional  "fitted  shock"  computer  code  was  develped  for 
the  first  time,  which  was  capable  of  predicting  a  flow  with  a  single  shock 
surface.  The  program  also  tracked  the  singularity  that  occurs  at  the  point 
where  the  impingement  shock  intersects  the  plume  boundary.  The  floating  shock 
program  became  far  to  complex  to  program  for  the  case  where  there  are  more 
than  a  single  shock  and  several  triple  points  in  the  flow.  A  different 
approach  was  investigated  which  does  not  explicitly  track  the  di scont i nui tes . 
The  idea  was  demonstrated  on  one  dimensional  unsteady  flow  to  trace  entropy 
discontinuities.  This  method  in  a  simple  way  tracks  a  slip  surface  using  a 
single  grid  point,  where  a  "capturing  method"  smears  the  discontinuity  over 
10-15  grid  points.  The  extension  of  these  ideas  to  three  dimensional  flows  is 
a  project  for  the  future. 

Theoretical  development  of  the  quantitative  description  of  three 
dimensional  flows  with  i nteresection  shock  surfaces  was  achieved.  The 
intersection  of  two  three  shock  surfaces  leads  to  a  complex  process.  The 
resulting  shock  pattern  depends  not  only  on  the  strength  of  the  two  shock 
waves  but  also  on  their  relative  orientation.  The  study  showed  how  a  local 
analysis  at  the  shock  intersection  line  could  be  used  to  explain  the  shock 
transmission/  reflection  configuration  at  the  intersection  line.  The  image  of 
this  intersection  process  in  the  hodograph  plane  (pressure/deflection  ) 
explains  the  requirement  for  transition  from  a  regular  reflection  process  to  a 
Mach  reflection  process. 
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The  flow  field  associated  with  the  axi  symmetric  plume  with  a  Mach  disc 
has  been  analyzed  using  a  new  iteration  scheme.  The  heart  of  this  approach  is 
the  integration  of  the  governing  interaction  equation  starting  from  the  sonic 
station  and  proceding  upstream.  This  analysis  which  leads  to  the  underlying 
interaction  equation  provides  theoretical  insight  for  the  first  time  into  the 
mechanisms  governing  the  Mach  disc  location.  The  analysis  also  includes 
viscous  mixing  effects  which  have  not  as  yet  been  studied  in  oast 
investigations.  The  integration  of  this  program  with  the  SPF  plume  code  is  a 
project  that  should  be  undertaken  in  the  future. 

Future  capabilities  in  prediction  of  three  dimensional  supersonic  flows 
as  complex  as  multinozzle  exhaust  will  increase  as  both  new  algorithms  are 
developed  and  computer  capabilites  increase.  Algorithm  development  must  be 
pursued  to  reduce  the  enormous  computer  logic  necessary  to  track  complex 
intersecting  shock  patterns.  Significant  computer  speed  and  size  developments 
could  relax  the  requirements  and  alter  the  shape  of  the  new  algorithms.  There 
is  a  great  amount  of  work  necessary  in  the  fundamental  understanding  of  three 
dimensional  flows.  There  are  a  variety  of  complex  conical  processes  such  as 
the  the  transition  of  the  shock  intersection  from  regular  to  irregular  and  the 
shock  lift  off  problem  that  must  be  studied,  understood  and  cataloged  in  order 
to  make  progress  and  achieve  prediction  capabilities  in  complex  three 
dimensional  supersonic  flows.  These  unit  solutions  are  the  three  dimensional 
counter  parts  of  the  familiar  two  dimensional  wedge  shock  solution  and  two 
dimensional  Prandtl  Meyer  solutions  that  are  used  extensively  as  building 
blocks  and  initial  conditions  in  two  dimensional  problems. 
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